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A.  BACKGROUND  -  PASSIVE  TARGET  TRACKING 


Passive  tracking  is  based  on  reception  and  analysis  of  signals  emitted  from 
a  target  at  an  unknown  time  and  location.  Range  measurement  to  the  target, 
usually  based  on  knowledge  of  the  travel  time  and  speed  of  the  signal,  is  not 
possible.  The  lack  of  range  information  makes  the  tracking  of  targets,  which  axe 
free  to  change  their  range,  a  challenge.  Our  interest  in  passive  target  tracking  is 
primarily  for  acoustic  tracking  of  targets  in  the  ocean.  However,  the  problem  of 
passive  localization  is  encountered  in  fields  like  radio  astronomy  and  seismology  as 
well  as  in  passive  sonar.  Several  passive  tracking  techniques  have  evolved  over  the 
years;  some  of  these  will  be  briefly  reviewed  here. 

1.  T.M.A  -  Target  Motion  Analysis 

Target  motion  analysis  is  a  range  estimation  technique  devised  for  World 
War  II  submarines.  Typical  targets  were  surface  ships  and  the  measurements 
were  primarily  sonar  bearing.  Occasional  course  and  speed  estimates  derived  from 
periscope  peeks  were  used  as  well,  but  mainly  for  initialization.  The  method  is 
based  on  hypothesizing  constant  course  and  speed  targets.  The  range  of  the  tar¬ 
get  whose  trajectory  would  have  produced  bearing  measurements  which  best  fit 
the  actually  observed  bearing  data  is  selected  as  the  estimate.  The  observability 
problem  is  partially  solved  by  maneuvering  the  observing  platform  while  assuming 
the  target  maintains  its  course  and  speed.  Own  ship  maneuvering  is  a  lengthy 
operation  and  the  assumption  is  that  the  target  maintains  course  and  speed  only 
of  limited  validity.  The  important  points  to  note  about  TMA  are: 


•  No  direct  or  indirect  measurement  of  range. 

•  Range  estimates  limited  to  nonmaneuvering  targets. 

•  Lengthy  own  ship  maneuver  requirement. 

•  Lengthy  solution  development  time. 

•  Very  accurate  bearing  measurements  needed. 

2.  Wavefront  Curvature 

Wavefront  curvature  is  an  indirect  range  measurement  technique.  It  is 
based  on  the  assumption  that  the  acoustic  waves  propagate  in  spherical  wavefronts. 
Signals  from  an  array  of  at  least  three  spatially  separated  receivers  are  delayed  in 
time  to  affect  focusing  .  The  amount  of  delay,  which  is  the  travel  time  difference 
of  arrival  (TDOA)  required  for  a  given  array  length,  is  then  translated  into  target 
range.  The  performance  depends  on  the  length  of  the  array  baseline,  the  precision 
of  both  the  sensor  location  along  the  array  and  the  time  delay  measurements. 

A  number  of  implementations  exist  with  array  baselines  varying  from  a 
few  tens  of  meters  in  integrated  systems,  to  hundreds  of  meters  in  large  multi¬ 
platform  distributed  configurations.  The  range  measurement,  which  is  indepen¬ 
dent,  can  however  be  further  combined  with  other  measurements  like  target  bear¬ 
ing  or  Doppler  in  order  to  reduce  the  ranging  error.  The  important  features  of  this 
method  are: 


•  Truly  independent  measurement  of  range. 

•  Requirement  fra  multiple  sensors  precisely  located  over  a  large  baseline. 


•  Applicability  to  maneuvering  targets. 

3.  Multipath  Tracking 


Multipath  (MP)  tracking  is  similar  in  principle  to  wavefront  curvature. 
Here  travel  time  differences  to  a  single  receiver  via  different  paths,  are  measured 


using  autocorrelation  techniques.  The  paths  most  commonly  used  are  the  direct 
path  and  those  reflected  from  the  surface  and  bottom  of  the  ocean.  The  reflected 
paths  can  be  viewed  as  paths  to  imaginary  virtual  receivers  in  mirror  image  lo¬ 
cations,  one  above  the  water  and  the  other  below  the  bottom.  In  this  sense,  the 
single  receiver  acts  like  an  array  of  three  receivers  (see  Section  C.2). 

As  in  the  wavefront  methods,  the  travel  time  differences  depend  on  the 
target’s  depth  and  range  and  on  the  array  baseline,  here  determined  by  the  water 
column.  If  propagation  is  assumed  to  be  along  straight  lines,  the  inverse  dependence 
of  target  position  on  travel  time  differences  is  relatively  simple.  It  provides  a 
mapping  of  time  delays  to  range  and  depth.  The  range  and  depth  observations 
thus  obtained  can  again  be  combined  with  other  available  target  measurements  in 
order  to  reduce  the  error.  The  main  features  of  multipath  tracking  are  therefore: 

•  Independent  measurement  of  range  and  depth. 

•  Compact  single  receiver  configuration. 

•  Applicability  to  maneuvering  targets. 

•  Dependence  on  ocean  propagation  conditions. 

The  MP  method  has  enjoyed  a  growing  amount  of  attention  in  recent  years 
due  to  the  advantages  stated  above  and  the  continued  improvement  in  time  delay 
estimation.  Compensating  for  the  dependence  of  the  method  on  ocean  conditions 
is  the  main  subject  of  this  work. 

B.  THE  BASIC  MULTIPATH  METHOD 
1.  The  Travel  Time  Differences 

The  acoustic  medium  of  the  ocean  is  confined  between  the  distinct  surface 
and  bottom  boundaries.  The  medium  thus  provides  multiple  reflective  paths  along 
which  sound  can  travel  between  two  given  points.  The  sound  travel  time  along 
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those  different  paths  will  differ  by  an  amount  dependent  in  a  known  manner  on 
the  geometry  of  the  scenario. 

Consider  the  direct  and  single  reflection  paths  shown  in  Fig.  1.1,  under 
the  following  assumptions: 


Fig.  1.1  Multipath  propagation. 

A.  The  sound  propagates  along  straight  lines. 

B.  The  bottom  and  surface  are  ideal  rigid  and  pressure  release  boundaries,  re¬ 
spectively. 

C.  The  travel  time  differences  can  be  analyzed,  resolved,  and  associated  with  the 
corresponding  paths. 

The  following  relations  then  exist  between  the  target  depth  D  and  range  R  position, 
measured  relative  to  the  observer,  and  the  travel  time  differences  rj  and  ti\ 
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(1.1a) 


n  =  T,  -  To  =  £  •  [(R2  +  D2  +  4 D2  -  4£>oD)1/2  -  P 
r2  =  r2  -  To  =  ^  [(R2  +  D2  +4(PU,  -  Do)2 -4(DW  -  D0)D)1/2  -  p](l.l&) 


Where  To,Ti,T2  axe  the  travel  times  of  the  direct,  surface,  and  bottom  paths, 
respectively;  Dw  is  the  water  depth  (  surface  to  bottom);  Da  is  the  observer  depth 
(surface  to  observer);  p  is  the  slant  range  given  by  (R2  +  D2)1/2  and  C  is  the  speed 
of  sound.  These  relation  can  be  inverted  to  express  the  slant  range,  the  depth  and 
the  range  as  functions  of  the  travel  time  differences.  The  results  are 

4  [(£>„  -  Do)  •  Dl  +  D,(DW  -  Z>„)!]  -  C2I?D„  +  T*(D„  -  D,)„  „  , 


2C[T1(D„-D„)  +  T!D„]  '  ' 

D = 4'(d;-po)  '  ^ + 2CT*P ~ 4(z>"  *  Do)^  (1-26) 


R=  (p2  —  D2)1*2 


(1.2c) 


Eq.  (1.1)  and  Eq.  (1.2),  referred  to  in  this  work  as  the  direct  and  inverse  functions 
respectively,  were  derived  by  Hassab  in  Ref.  1. 

2.  Time  Delay  Estimation 


The  signal  y(t),  which  is  the  signal  y'(t)  time  shifted  by  the  travel  time  of  the  first 
arrival  To,  that  is, 

y(t)=y'(t-T0)  (1.4) 


can  be  expressed  in  terms  of  the  TDOAs  r,-  as 

N  N 

y(t)  =  ~T°  ~Ti)-  -  r«)  (1.5) 

i=0  1=0 


where  t*  =  T{  —  To.  If  the  source  of  the  target’s  acoustic  emission  is  propeller 
cavitation  noise,  as  is  the  case  for  broadband  MP,  then  the  original  signal  x(t)  is 
a  broiwuuond  random  process.  The  autocorrelation  of  its  received  version  y(t)  will 
exhibit  peaks  at  time  lags  equal  to  r,.  The  autocorrelation  can  be  easily  computed 
at  discrete  time  lags  using  a  sampled  version  of  the  received  signal.  As  long  as 
the  ACF  peaks  are  resolvable,  they  can  be  interpolated  between  the  discrete  lags 
providing  the  required  continuous  time  delay  measurement.  The  measured  delays 
can  be  further  smoothed  using  a  time  delay  tracker. 

In  the  narrowband  case,  low  frequency  tonals  emitted  by  target  machinery 
can  be  received  at  long  range.  Since  the  autocorrelation  function  of  such  sinusoids  is 
periodic,  the  ACF  approach  is  impractical.  A  different  approach  has  been  applied 
to  this  situation  which  is  based  on  standing  waves  generated  by  the  tonals  in 
the  water  column.  The  sound  field  intensity  is  sensed  by  a  number  of  receivers 
placed  along  a  vertical  array.  A  technique  called  field  matching  is  used  to  select  a 
target  position  which  will  give  rise  to  a  sound  field  that  best  matches  the  measured 
intensity.  The  effort  in  this  work  concentrates  on  the  broadband  case. 


The  MP  depth  and  range  measurements  are  contaminated  with  noise. 
Sources  of  the  noise  include  acoustical  noise  which  distorts  the  autocorrelation 


of  the  received  signal,  errors  and  fluctuations  of  the  MP  conditions  in  the  ocean, 
and  errors  in  the  receiver  and  signal  processor.  To  estimate  the  actual  position 
from  the  noisy  measurements,  a  target  tracker  is  employed.  The  tracker  fits  the 
measurements  with  an  estimate  that,  on  the  average,  minimizes  the  squared  error 
between  the  actual  and  estimated  position.  Such  a  tracker,  also  referred  to  as 
an  estimator,  observer,  or  filter ,  may  also  use  other  available  measurements  like 
bearing  and  Doppler. 

Many  trackers  of  the  type  described  above  are  available.  These  range 
from  simple  averaging  trackers  to  large  banks  of  Kalman  filters.  The  trackers  differ 
primarily  in  the  measurements  they  use  and  the  complexity  of  the  assumed  target 
model. 


A  typical  complete  MP  measurement  and  tracking  system  is  shown  in 
Fig.  1.2.  Estimated  TDOAs  are  converted  to  depth  and  range  by  the  prefilter. 
Depth  and  range  are  then  combined  with  bearing  to  form  a  3-D  target  position 
measurement  which  is  filtered  by  the  target  tracker. 

Current  MP  tracking  methodology  is  lacking  in  the  following  areas: 

•  The  methods  assume  straight  line  propagation  which  is  true  only  in  a  homoge¬ 
neous  medium.  This  assumption  is  largely  in  error  for  the  realistic  inhomoge¬ 
neous  (IH)  case.  This  causes  large  errors  in  the  transformation  of  time  delays 
to  depth  and  range. 

•  The  assumption  that  path  and  delays  are  easily  associable  is  wrong  for  many 
practical  cases.  A  significant  amount  of  detailed  oceanographic  data  is  re¬ 
quired  in  order  to  identify  the  multiple  paths.  This  information  is  not  typically 
available  in  an  actual  target  tracking  situation. 

•  Receiver  time  delay  resolution  is  not  perfect  but  limited  by  the  bandwidth 
of  the  signal  and  the  receiver.  This  affects  the  performance  whenever  the 
geometry  yields  two  or  more  multipaths  with  similar  travel  times.  Three  such 
typical  geometries  Eire  : 
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Fig.  1.2.  MP  tracking  system, 
o  Long  range,  where  all  paths  tend  to  be  of  same  length. 

o  Target  or  observer  close  to  the  surface  or  bottom  where  the  direct  path 
and  the  reflected  path  have  similar  length. 

o  Target  and  observer  at  symmetrically  opposing  depths  such  that  the  bot¬ 
tom  and  surface  paths  have  similar  travel  times. 

•  The  lack  of  exact  knowledge  of  bottom  depth  and  structure  limits  both  the 
time  delay  measurement  and  the  exact  translation  to  depth  and  range. 

•  Other  effects  including  slanted  bottom,  and  the  anisotropic  nature  of  the  bot¬ 
tom  which  may  have  different  slopes  at  different  directions. 

Our  work  seeks  to  address  some  of  these  effects,  in  particular  the  first 
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In  1976  Hassab  published  an  original  work  titled  “Passive  tracking  of 
moving  source  by  a  single  observer  in  shallow  water”  [Ref.  1]  which  was  the  starting 
point  for  the  research  in  the  area.  In  this  work,  the  use  of  measured  travel  time 
differences  as  inputs  to  a  target  tracking  filter  was  established.  The  tracker  was 
realized  in  cartesian  coordinates  using  an  extended  Kalman  filter  and  assuming  a 
nonmaneuvering  target.  Earlier  work  by  Hassab  [Ref.2]  dealt  with  the  very  process 
of  measuring  the  travel  time  differences  using  a  technique  based  on  the  cepstrum. 

b.  Singer 

In  the  early  T.M.A  algorithm  the  development  of  a  range  solution  was 
halted  when  a  target  maneuver  was  detected  (a  “Zig  Zag”  in  World  Welt  II  jargon). 
The  maneuvers  were  detected  when  constant  tracking  errors  indicated  a  mismatch 
between  the  actual  target  and  its  model.  With  the  application  of  Kalman  filters 
to  target  tracking,  the  maneuvering  command  uncertainty  was  represented  as  a 
white  gaussian  process  noise.  Singer  [3]  improved  the  model  by  coloring  (low  pass 
filtering)  the  command  noise  to  correspond  to  the  expected  maneuver  dynamics. 

c.  Moose 

In  the  1970’s  multiple  model  (multiple  hypothesis)  methods  were  in¬ 
troduced  [4]  to  address  the  maneuvering  problem.  Here  a  number  of  target  models 
based  on  different  hypothesized  maneuvers  are  computed  in  parallel,  and  combined 
in  a  Bayesian  manner  to  form  the  overall  estimate. 

This  approach  was  applied  to  both  airborne  and  underwater  targets  by 
Moose  [5,11]  and  his  colleagues  McCabe  [6],  Gholson  [7],  Van  Landingham  [8],  Dai¬ 
ley  [9],  and  Caputi  [10].  Three  main  issues  related  to  MP  tracking  were  addressed 


in  this  work,  which  is  ongoing;  They  are  the  selection  of  a  particular  coordinate 
system,  the  account  for  target  maneuver,  and  the  bias  associated  with  the  nonlin¬ 
ear  transformation  of  TDOAs  to  depth  and  range.  Performance  evaluations  were 
carried  out  in  Moose’s  work  but  both  the  simulation  and  tracking  were  done  with¬ 
out  taking  account  for  the  effects  of  the  IH  acoustic  medium  and  the  realistic  delay 
estimation  process. 

2.  Time  Delay  Estimation  and  Source  Localization. 

The  area  of  time  delay  estimation  has  received  enormous  attention  in  the 
last  twenty  years.  Among  the  contributors  are  Schultheiss  [13],  Ianniello  [14],  Wien- 
stien  [15],  FYiedlander  [16,  17],  and  others.  Work  in  this  area  covered  the  topics 
of  estimation  instrumentation,  theoretical  and  experimental  error  bounds,  estima¬ 
tion  resolution,  continuous  delay  reconstruction  from  sampled  time  sequences,  and 
more.  A  very  good  summary  is  presented  in  Refs.  15  and  16. 

Use  of  TDOA  for  source  localization  was  investigated  by  Schultheiss  for  a 
variety  of  receiving  array  configurations.  In  his  work  [18]  lower  bounds  on  achiev¬ 
able  error  using  bandlimited  signals  is  developed  using  the  Cramer  Rao  lower  bound 
(CRLB)  and  Ziv  Zakai  lower  bound  (ZZLB).  In  the  Ph.D.  dissertation  by  his  stu¬ 
dent  Hamilton  [19]  the  lower  bound  for  the  MP  localization  error  variance  is  also 
developed.  Hamilton  also  establishes  the  relation  between  the  wavefront  curvature 
and  the  multipath  ranging  using  reflected  mirror  images  of  the  receiver. 

3.  Acoustic  Propagation  and  Modeling 

The  basic  form  of  sound  propagation  in  an  inhomogeneous  medium  where 
the  speed  of  sound  is  a  function  of  depth  has  long  been  known.  The  development 
of  the  associated  approximations  and  the  resulting  ray  acoustics  are  presented,  for 
example,  in  Ref.  20  by  Ziomek. 
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Numerous  computer  models  have  been  developed  to  compute  sound  prop¬ 
agation  in  the  ocean  based  on  ray  acoustics,  starting  from  analog  computer  sim¬ 
ulations  which  were  used  for  ray  tracing,  and  proceeding  up  to  very  large  and 
sophisticated  numerical  models  which  account  for  many  of  the  special  effects  of 
the  ocean  medium.  A  classic  model  in  this  category  is  the  Generic  Sonar  Model 
developed  by  Weinberg  [21]. 

Only  a  few  models  address  the  problem  of  finding  the  rays  traveling  be¬ 
tween  two  given  end  points,  the  eigenray  problem.  A  recent  model  of  this  kind  is 
SMART  (SMall  Acoustic  Ray  Tracer)  developed  by  Novick  [22]. 

Research  on  the  nature  of  reflection  from  the  bottom  and  the  surface  of 
the  ocean  also  dates  back  to  World  War  II.  Contributors  in  this  field  include,  for 
example,  Clay  and  Medwin  [23].  Recent  work  in  this  area  reveals  the  interrela¬ 
tionship  between  the  various  oceanographic  processes,  for  example,  internal  waves, 
surface  interaction  of  the  ocean  and  the  atmosphere,  and  shear  waves  induced  in 
the  bottom  of  the  ocean  by  a  sound  field.  The  limitation  of  simplified  partial, 
lumped  models  for  these  phenomena  is  becoming  apparent. 

D.  PROBLEM  STATEMENT 

This  thesis  research  seeks  to  address  the  following  issues: 

•  To  account  for  the  complex  effects  of  the  inhomogeneous  (vertically  stratified) 
ocean  medium  on  broadband  MP  tracking. 

•  To  investigate  the  impact  and  account  for  the  effects  of  a  realistic  receiver  and 
time  delay  estimator  on  the  overall  tracking  process. 

In  addition  our  work  includes: 

•  Modification  of  the  state  of  the  art  target  tracking  algorithm  to  decrease  some 
of  its  estimation  bias. 


•  Evaluation  of  the  complete  tracking  algorithm  using  a  realistic  environment, 
and  target  simulation  considering  all  main  sources  of  estimation  errors. 

The  main  contributions  and  key  discussions  of  each  of  the  foregoing  issues  are  given 

below. 

1.  The  Inhomogeneous  Ocean  Medium 

The  real  ocean  is  acoustically  inhomogeneous  ,  in  that  among  other  effects, 
the  speed  of  sound  varies  with  depth.  This  variation  significantly  effects  the  MP 
propagation  and  travel  time  as  shown  in  Fig.  1.3  where  the  receiver  ( Rx )  is  at 
range  =  0.  Note  that  the  direct  path  between  the  source  and  the  receiver  is  com¬ 
pletely  eliminated  due  to  the  ray  bending  in  Fig.  1.3b  (am  exhaustive  search  for  the 
eigen  rays  was  conducted  for  this  plot  and  the  resultant  eigenrays  are  plotted).  The 
effects,  believed  to  be  analyzed  quantitatively  here  for  the  first  time,  are  shown  to 
render  the  assumption  of  straight  line  propagation  to  be  of  limited  practical  use. 
Accounting  for  this  effect  is  difficult  since  the  inverse  function,  which  transforms 
travel  time  differences  to  depth  and  range  is  not  readily  computable  by  existing 
acoustic  models.  Devising  an  inversion  method  to  account  for  the  IE  effects  on  the 
MP  tracking  is  a  main  goal  of  this  research. 

The  acoustic  energy  attenuates  as  it  propagates  through  the  water,  and 
many  factors  contribute  to  this  loss  including  spreading,  absorption,  and  reflec¬ 
tion.  While  these  loses  and  their  impact  on  the  time  delay  estimation  noise  are 
relatively  well  understood,  they  were  nevertheless  never  considered  in  past  simula¬ 
tions.  Inclusion  of  the  reflection  and  spreading  effects  to  increase  the  reliability  of 
the  simulation  is  another  part  of  the  effort  to  better  account  for  the  effects  of  the 
medium. 
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Fig.  1.3.  IH  MP  propagation,  a.  with  direct  path,  b.  without  direct  path 


2.  Realistic  Receiver  and  Delay  Estimation 

a.  Limited  Resolution 

The  finite  bandwidth  of  the  target  signal  and  receiver  limit  the  time 
delay  resolution  to  a  few  tenths  of  a  millisecond.  A  significant  amount  of  this 
research  was  dedicated  to  prediction  and  improvement  of  this  limitation  but  less 
work  was  aimed  at  investigating  the  details  of  its  impact  on  MP  tracking.  The 
effect  of  this  limited  resolution  on  the  MP  tracker  and  the  means  to  overcome  it 
are  another  subject  addressed  by  this  research. 

b.  Nonidentiflable  Paths 

The  polarity  of  the  ACF  can  support  association  of  the  time  delays 
with  their  corresponding  paths  when  the  MP  structure  is  known  to  be  simple. 
However,  when  the  structure  is  complex  and/or  unknown,  such  as  in  cases  where 
there  is  a  lack  of  direct  path  due  to  ray  bending,  the  simple  association  is  not 
possible.  The  ability  to  associate  delay  with  paths,  assumed  in  previous  work,  is 
not  assumed  here. 

c.  Nonnegative  Noise 

The  even  symmetry  of  the  autocorrelation  function  does  not  enable 
the  measurement  of  time  delay  polarity.  In  the  multipath  situation,  there  is  a  first 
arrival  followed  by  lagging  replicas.  This  justifies  the  use  of  the  positive  time  delays 
only.  The  implication  is  that  the  delay  noise  is  not  normally  distributed  as  was 
assumed  in  the  past.  A  distribution  which  is  zero  for  negative  delay  is  a  better 
description  and  it  leads  to  a  nonzero  delay  estimation  bias  error. 

3.  Tracker  Modification 

In  the  maneuvering  target  tracker  introduced  by  Moose,  the  conditional 
mean  of  some  hypotheses  commands  is  used  as  the  estimate.  This  estimate  is 
biased  if  the  hypotheses  commands  are  not  symmetrically  centered  around  the 
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actual  value.  Removing  this  bias  is  difficult  since  modification  of  the  hypotheses 
disrupts  a  recursion  used  in  the  estimation  process.  Modification  of  the  tracker  to 
remove  the  bias  without  disrupting  the  recursion  or  introducing  large  estimation 
transients  is  one  of  the  goals  of  the  research. 

4.  Evaluation 

Many  elements  are  involved  in  MP  tracking,  namely 

•  The  target  and  its  dynamics. 

•  The  medium  and  its  multipath  structure. 

•  The  receiver  and  the  time  delay  analyzer. 

•  The  conversion  from  time  delays  to  position  measurement. 

•  The  actual  acquisition  and  tracking  of  a  maneuvering  target. 

Evaluation  of  the  system  performance  is  done  in  two  steps.  First,  each  and 
every  component  is  evaluated  separately.  This  is  done  under  simulated  conditions, 
which  are  determined  for  each  component  based  on  the  analysis  of  the  overall 
system  operation.  Then  and  only  then,  the  integrated  system  is  evaluated,  using 
the  results  of  the  individual  component  evaluation  as  ‘reference  data’. 

E.  SCOPE  AND  OUTLINE 

This  work  approaches  the  MP  tracking  as  a  system  problem.  The  remaining 
chapters  are  as  follows. 

Chapter  Two  deals  with  the  3-D  target  tracker.  A  state-of-the-art  tracker 
is  described.  An  idealized  receiver  delay  analyzer  and  homogeneous  straight  line 
propagation  axe  assumed.  The  algorithm  provides  a  good  starting  point  in  terms 
of  its  treatment  of  a  3-D  maneuvering  target.  New  interpretation  and  analysis  is 
provided  for  certain  aspects  of  the  estimation  mechanism,  leading  to  a  modification 
and  improvement  .  A  systems  approach,  not  used  in  previous  work,  is  applied  to 


evaluate  the  performance  of  the  modified  tracker.  Some  of  the  parameters  and 
algorithms  used  for  the  simulation  are  not  discussed  in  detail  here,  since  they 
result  from  the  removal  of  the  above  idealized  assumptions  which  is  explained  in 
Chapter  Three. 

Chapter  Three  addresses  the  main  source  of  tracking  error,  namely  the  IH 
medium  together  with  the  effects  of  a  realistic  receiver  and  delay  analyzer.  The 
nature  of  the  problem  is  described  and  the  proposed  solution  is  discussed  in  detail. 
The  performance  of  the  method  is  then  evaluated  and  analyzed  in  detail. 

Chapter  Four  presents  a  detailed  evaluation  and  analysis  of  the  performance 
of  the  new  inversion  prefilter. 

Chapter  Five  resumes  the  overall  system  performance  analysis,  this  time  with¬ 
out  idealized  assumptions  and  using  the  new  inversion  algorithm  and  the  improved 
tracker. 

The  work  is  summarized  in  Chapter  Six  where  conclusions  and  recommended 


extensions  are  presented. 


A.  INTRODUCTION 


The  three  dimensional  (3-D)  target  tracker  developed  in  this  chapter  is  based 
on  a  state-of-the-art  design  by  Moose,  and  McCabe.  The  modified  Kalman  filter 
which  forms  the  central  portion  of  the  tracker  was  recently  investigated  by  Saez 
[Ref.  24],  in  conjunction  with  this  study.  In  his  work,  Saez  includes  a  detailed 
review  of  the  tracker,  results  of  which  are  briefly  repeated  here  for  completeness 
and  to  establish  notation. 

The  emphasis  in  this  section  is  on  improvements  of  the  tracker  design  and  a 
new  overall  systems  approach  to  the  performance  evaluation.  The  modifications  in¬ 
clude  an  adaptive,  instead  of  a  fixed,  command  hypothesis  bank  and  an  advancing 
smoother.  Both  modifications  are  intended  to  reduce  tracking  biases.  The  investi¬ 
gation  also  includes  the  use  of  a  second  instead  of  a  third  order  target  model.  This 
was  found  to  reduce  computational  load  without  sacrificing  estimation  accuracy. 

A  realistic  new  model  for  the  time  delay  estimation  noise  is  used  in  evaluating 
the  performance.  The  model  incorporates  propagation  effects  as  well  the  effects 
of  some  inaccuracies  in  the  time  delay  estimation.  This  enables  a  more  realistic 
evaluation  of  the  overall  performance  for  the  homogeneous  case,  which  will  be 
extended  in  Chapters  Three  through  Five  to  the  IH  case. 

A  range  and  be  airing  coordinate  decoupling  approximation  was  introduced  in 
the  original  tracker  by  McCabe  [Ref.  6]  to  linearize  the  model.  An  interesting 
interpretation  of  this  procedure  provides  an  explanation  for  errors  in  tracking  non¬ 
maneuvering  targets  that  occur  at  short  ranges. 


B.  TARGET  DESCRIPTION 
1.  Horizontal  Plane 

A  second  order  physical  model  is  used  to  describe  target  motion  in  the  X 
and  Y  directions  of  the  horizontal  plane.  Controlling  the  motion  is  a  command 
which  forms  the  third  degree  of  freedom  in  each  axis.  If  Th(t )  is  the  command 
thrust  in  Newtons  (N)  and  Dc  is  the  drag  in  [N  per  m/sec]  the  differential  equation 
for  each  axis  which  follows  from  Newton’s  second  law  is 


Th{t )  —  Dcx(t)  =  mx(t).  (2.1) 

This  equation  can  be  rewritten  sis 

x(t)  =  —  a(U(t)  —  x(<))  (2.2) 

where  a  =  Dc/m  [sec-1]  is  the  reciprocal  of  the  system’s  damping  time  constant, 
and  U  =  Th(t)/Dc  [m/sec]  is  the  speed  command,  i.e.,  the  speed  at  which  the  plat¬ 
form  will  move  when  steady  state  is  reached.  If  it  is  assumed  that  X(0)  =  X(0)  =  0 
and  that  the  command  is  a  constant  U  then  the  solution  is  : 

i(t)  =  U(l  -  e'at)  (2.3) 
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2.  Depth 


Because  of  the  reduced  dynamics  expected  for  underwater  targets  in  the 
depth  channel,  the  target’s  motion  in  depth  was  modeled  as  a  first  order  system. 
The  command  Ud  [m]  can  therefore  be  represented  in  terms  of  the  steady  state 
depth,  and  the  differential  equation  along  the  depth  axis  is 

D(t)  =  ad(Ud-D(t)).  (2.5) 
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The  solution  of  Eq.(2.5)  is  given  by 


D(t)  =  D{o)e~ait  +  Ud  •  (1  -  e-Q-‘) .  (2.6) 

3.  Command  Noise  Augmentation 

Eq.  (2.2)  describes  target  motion  in  response  to  a  known  control  U .  How¬ 
ever  the  control  (command)  is  not  known  to  the  observer,  especially  if  the  target  is 
maneuvering.  Singer  partially  accounts  for  this  by  adding  a  colored  noise  compo¬ 
nent  w'  to  the  command.  This  noise  w'  is  modeled  by  a  first  order  recursive  filter 
driven  by  a  normally  distributed  white  noise  input  w. 

=  -awxv'(t)  -I-  w(t)  (2.7) 

The  lowpass  model  was  chosen  since  it  represents  a  maneuver  which  typically  takes 
at  least  a  few  seconds  to  complete.  The  command  is  thus  modeled  in  the  ob¬ 
server  by  the  sum  of  an  assumed  command  U  and  the  random  noise  w' ,  that  is, 
Uobserver  —  U  +  w'(t).  This  description  suggests  estimation  of  the  command  U  by 
the  multi-model  (MM)  estimator  described  in  Section  C. 

In  order  to  maintain  a  complete  system  state  description,  the  coloring 
LP  filter  and  its  state  w'  are  combined  into  the  target  state  equations.  Both  the 
command  U  and  the  the  noise  w  are  represented  in  units  of  the  resulting  steady 
state  speed.  The  differential  equations  are  then  represented  in  discrete  form  by  the 
following  set  of  (3-D)  Cartesian  motion  difference  equations. 


X„  =  ^X„_i  +  TUx  +  n— 1 
Yn  =  Y„^Y„_i  +  TUy  +  Vtwy 
Dn  =  ad  Dn~ i  +  (1  —  ad)Ud 
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(2.8a) 

(2.86) 

(2.8c) 


where 


X  =  (z  x  w'z)T 
Y  =  (y  y  w',)T 

and  the  components  of  the  matrices  I\  ¥  are  detailed  in  Appendix  A,  and 
ad  =  e~°iT.  Note  that  since  the  equations  are  linear  and  decoupled,  a  change 
in  target  course  at  a  constant  speed  is  modeled  here  as  an  acceleration  along  one 
axis  and  a  deceleration  along  another. 

4.  Cylindrical  Coordinate  Observations 

Most  surveillance  and  tracking  sensors  produce  measurements  relative  to 
their  own  position  and  orientation.  The  same  is  true  in  the  MP  case,  where  range 
and  depth  axe  indirectly  deduced  from  the  travel  time  differences  across  the  vertical 
plane.  Bearing  is  measured  by  beamforming;  here  sonar  azimuthal  beams  are 
formed  using  travel  time  differences  across  the  horizontal  plane.  Thus,  the  natural 
measurements  are  in  the  cylindrical  coordinate  system  of  range,  depth  and  bearing. 
These  are  relative  to  the  observing  sensor  wnich  is  therefore  set  at  the  origin.  The 
detailed  set  of  variables  used  to  describe  the  positions  of  the  observer  and  the  target 
in  the  ocean  is  shown  in  Fig.  2.1  and  defined  in  Table  2.1. 

TABLE  2.1 

CYLINDRICAL  COORDINATES  VARIABLES 

Do  Observer  depth  measured  from  surface. 

Dw  Depth  of  the  water. 

D  Depth  of  target  relative  to  observer. 

R  Horizontal  range  of  target  relative  to  observer. 

B  Bearing  (azimuthal  target  angle  relative  to  north). 

B  Bearing  rate. 

Ur  Speed  command  in  range  direction. 

Ub  Speed  command  in  cross  range  direction. 

U d  Depth  position  command. 

U'd  Depth  speed  ctd{Ud  —  D) 

p  Slant  range. 
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Fig.  2.1  Cylindrical  coordinates. 
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5.  Linearization  and  Decoupling 

In  order  to  use  the  cylindrical  measurements  as  observations  for  the  Carte¬ 
sian  target  model  of  Eq.  (2.8),  the  measurements  have  to  be  transformed.  This 
leads  to  a  coupled,  nonlinear  observer.  Two  possible  approaches  have  been  applied 
to  this  problem  in  the  past.  One  is  based  on  conversion  and  linearization  of  the 
measurements,  and  leads  to  the  extended  Kalman  filter,  (EKF).  This  approach  was 
applied  to  MP  tracking  by  Hassab  [1].  The  other  was  introduced  by  Moose  [5]  and 
is  based  on  a  decoupling  and  linearization  of  the  state  equations,  using  discretized 
approximations . 

Investigation  of  the  two  approaches  [7]  indicated  a  preference  for  the  second 
since  taking  account  of  target  maneuvers  in  the  EKF  seemed  to  require  a  very  large 
computational  load.  The  second  approach,  developed  in  detail  by  McCabe  in  [6], 
is  based  on  the  following  main  steps. 

The  Cartesian  target  state  description  is  transformed  to  cylindrical  coor¬ 
dinates  using  the  transformation: 


R  =  (x*  +  y*)1'2 

(2.9  a) 

B  —  tan-1  ( y/x ) 

(2.96) 

D  =  D. 

(2.9c) 

In  the  bearing  channel,  the  range  is  assumed  constant  for  the  duration  of 
the  sampling  period.  The  speed  maneuver  commands  Uz  and  Uy  along  the  X  and 
Y  directions  are  replaced  by  commands  along  the  range  and  cross  range  (bearing) 
coordinates  Ur  and  Ub  respectively. 

Range  and  bearing  coordinate  coupling  is  thus  reduced  to  a  parametric 
relation.  This  requires  the  bearing  channel  matrices  to  be  recomputed  only  when 


the  range  changes  by  more  then  some  given  ratio  and  no  more  then  once  per 
iteration.  The  resulting  state  equations  are 


Dn  =  adDn-i  +  (1  —  ad )  U dn_i 

R-n  =  ^r^n- 1  "b  T'rUr  n  — 1  "b  ^ RWr 

Bn  =  ^Bn_,  +  TbUb  n  — X  +  VbWb 


(2.10a) 

(2.106) 

(2.10c) 


where 


R  =  {R  R  w'r)T 

(2.11) 

h 

and 

«\N 

B  =  (B  B  w'„)T. 

(2.12) 

it. 

Each  channel  of  the  three  decoupled  channels  (depth,  range,  bearing) 

is  observed 

g 

by  a  scalar  z  defined  as 

e.  ^ 

cV 

• 

;  •  * 

Zdn  —  d„  -f*  t’o|n 

(2.13a) 

>: 

Zm  —  Rn  "h 

(2.136) 

v- 

*-■> 

%bn  —  Bn  -f-  Vbn 

(2.13c) 

s? 

■~vi 

where  vd ,  vr,  and  vb  are  the  observation  noises.  The  components  <f>b  Tj  9b  of 

•V« 

.*>y 

the  cross  range  matrices  are  range  dependent  and  presented  in  Appendix  A.  The 
subtle  terminology  -  preference  of  the  term  “cross  range”  to  “bearing”  -  used  by 
McCabe,  is  significant.  It  is  indicative  of  the  fact  that  while  the  bearing  variable 
is  used,  a  constant  cross  range  rate  motion  is  actually  being  modeled  as  the  basic 
nonmaneuvering  case.  The  command  Ub  is  a  speed  command  in  the  cross  range 
direction  and  not  a  bearing  rate  command.  While  the  former  is  a  constant,  the 
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latter  is  range- dependent.  Similarly  the  range  rate  command  Ur  reflects  an  as¬ 
sumed  constant  range  rate  for  the  nonmaneuvering  case.  This  again  is  an  incorrect 
assumption  for  a  target  moving  along  any  straight  line  other  than  the  line  of  sight. 
The  impact  of  these  assumptions  will  become  more  apparent  in  the  next  section. 

C.  MULTI-MODEL  ESTIMATION 
1.  Concept 

Use  of  the  classic  Kalman  filter  as  an  observer  for  a  maneuvering  target 
is  possible  only  if  provision  is  made  for  the  unknown  maneuver  command.  The 
method  suggested  by  Singer  \ct  i  not  suffice  in  itself  since  it  dictates  maintaining 
high  Kalman  gain;  this  inhibits  effective  measurement  noise  filtration. 

A  better  procedure  is  known  as  the  Multi- Model  (MM)  [43]  and  was  first 
applied  to  the  MP  tracking  problem  by  Moose.  In  the  MM  a  number  (N)  of 
commands  are  hypothesized  to  form  a  command  bank  vector  UI  defined  as: 

V1  =  (UUU2.-.UN).  (2.14a) 

The  bank  is  ideally  centered  around  the  mean  command  value  and  the  commands 
in  the  bank  are  evenly  spread  to  span  the  full  range  of  possible  commands.  A 
corresponding  bank  of  Kalman  filters  (models)  is  formed  and  each  filter  in  the 
bank  separately  tracks  the  target  using  its  particular  hypothesized  command  and 
the  common  measurement  Z.  The  estimated  states  Xi  from  all  the  filters  in  the 
bank  are  combined  to  form  a  conditional  mem  which  is  used  as  the  overall  estima¬ 
tor  output.  The  conditional  mem  depends  on  the  discrete  commmd  probability 
distribution  represented  by  the  weight  vector  W  defined  as 

w  =  (W1  W2  ...Wn)t  (2.146) 
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where  the  weight  Wi  is  the  conditional  probability  (weight)  that  a  command  U,  was 
exercised  in  the  previous  time  period  given  all  past  observations  zn.  The  expression 
for  the  estimated  state  X  is  £,=i  XiWi. 

The  discrete  distribution  of  W  can  be  computed  recursively  if  the  com¬ 
mand  is  modeled  as  a  semi-markov  process  with  a  known  probability  transition 
matrix  9.  The  conditional  probability  of  each  command,  given  the  past  position 
measurements  and  the  hypothesis,  is  computable  using  the  Kalman  filter’s  error 
propagation  matrix  P.  A  diagonal  matrix  A  with  elements  a[i,  ij  set  proportional 
to  the  conditional  probability  of  the  innovation  given  the  hypothesis  i  is  also  used 
in  the  recursion  which  is  given  below.  An  overall  block  diagram  of  the  resulting 
estimator  is  shown  on  Fig.  2.2.  The  development  of  the  basic  estimator  is  reviewed 
in  detail  in  Ref.  24  and  will  not  be  repeated  here  except  for  the  resulting  estimator 
equations.  Emphasis  in  the  following  sections  will  be  directed  towards  the  more 
recent  evaluation  and  modification  of  the  algorithm  and  to  its  application  to  the 
IH  multipath  tracking  problem. 


2.  The  Multi-Model  Estimator  Equations 


Fig.  2.2  Multi-Model  Block  Diagram 
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TABLE  2.2 

MM  VARIABLE  DEFINITIONS 

3x3  state  transition  matrix. 

3x1  control  gain  vector. 

3x1  Singer  process  noise  gain  vector. 

3x1  observation  matrix  H  =  [1,0,0]. 

3x1  Kalman  gain  matrix. 

3x1  state  vector  X  =  ^  j  at  time  n. 

3x1  estimate  of  filter  t  in  the  bank  at  time  ni 
given  past  observations  z\ ,  z-j.  . . .  zn 2 
The  observation  scalar. 

Number  of  models  in  the  bank. 

1  x  N  a  vector  of  l’s. 

3  x  N  state  bank  matrix.  XI [t  j]  is  the  i-  th  component 
of  the  state  vector  X;  of  the  j-th  model, 
i.e.,  XI  =  [XijXjj  •  -|X/v].  The  matrix  is  time 
index  like  Xj  above. 

The  assumed  command  of  the  ith  filter  in  the  bank. 

1  x  N  command  bank,  evenly  spanning  the  Umm  —  Umaz 
interval  of  allowable  commands,  UI  =  [{/j ,  U? , . . .  U/v ] . 

Minimum  command  in  the  bank. 

Maximum  command  in  the  bank. 

Command  bank  separation  Umax  —  Umi„. 
lxl  command  bank  step  size  UI[2]  —  UI[1]. 

Center  of  command  bank  {Umaz  —  Um «n)/2. 

lxl  zero  mean  gaussian  measurement  noise  v  ~  N(0,crl). 

lxl  measurement  noise  standard  deviation. 

lxl  Singer  white  process  noise  input. 

lxl  standard  deviation  of  the  conditional  innovation 

distribution  given  the  hypothesis. 

Standard  deviation  of  the  command  quantization  error. 

3x3  command  quantization  error  matrix. 

N  x  1  hypothesis  command  probability  (weight)  vector. 

W  =  [Wx  W2  . . .  WN] 

N  x  N  command  Markov  probability  transition  matrix,  flu  jj 
-  the  probability  that  command  j  will  change  to  command  1  in  one  step. 

d[i,j]  =  (  1  f? r  -  \  *  ,  ^  1  where  Pr  is  the  probability  of 

l  TT^Tv1  -  yr)  1  T  J  J 
unchanging  command. 

N  x  N  conditional  innovation  probability  matrix. 

Smoothed  output  of  the  estimates  X,f7,  W 
Smoothing  coefficient 


The  estimation  equations  are  presented  in  the  order  of  their  use  in  the 


recursion.  The  Kalman  filter  equations  are 

XIn,n_1  =^.XIn_1|n.1+r.UI  (2.15) 

P»|n-1  =  +  Du  +  (2.16) 

Gn  =  Pnin-iH7,  [HPnln.1HT  +  *1]  (2.17) 

Pn|n  =  [I  -  G„H]  Pn|n_!  (2.18) 

XIn|n  =  XInjn_!  +  G„  (r„  •  N  —  H  •  XIn|n_ i )  (2-19) 


Note  that  the  state  of  the  multiple-model  XI  is  the  3  x  jV  state  bank  matrix  whose 
columns  represent  the  states  of  the  individual  filters  in  the  bank.  The  corresponding 
command  UI  is  an  N  x  1  row  vector  which  represents  the  entire  command  bank. 
D.  is  the  command  quantization  error  matrix  and  is  computed  once  as 

AU2  t 

D„=r— rr.  (2.20) 

Since  Eq.(2.15)  through  (2.19)  are  common  to  all  the  filters,  the  classical  single 
channel  Kalman  gain  vector  Gn  can  be  computed  once  and  used  for  all  the  filters. 
The  adaptive  command  estimation  equations  are 


^/n  —  HP ra|n_1H:r  -f-  a\ 

(2.21) 

fe  -( z„  N— H  XI(l,»])J/2<7/«  j  _  'I 

[,’;ln  “  l  0  «  7*  j  J 

(2.22) 

W'n  =  (N  A  0  Wn_,]_1  A0W „_i 

(2.23a) 

Wn  =  au,Wn_1  +(l-aw)W'n 

(2.236) 
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where  the  vector  of  conditional  command  probability  W  is  recursively  computed. 
The  conditional  command  probability  W  is  then  used  to  compute  the  conditional 
mean  estimates  X  and  U  of  the  state  and  command  banks  XI  and  UI  respectively. 


Un  =  UI  W„  (2.24) 

X„  =  XI„  •  Wn  (2.25) 

Additional  first  order  smoothing  is  applied  to  the  MM  estimate  in  the  form  of 

(Xop)n  =  Ur (X0p)n_1  +  (1  -  az)Xn  (2.26) 

( Uop)n  =  au(Uop)n- 1  +  (1  -  au)Un  (2.27) 

(Wop)n  =  att(Wop)n_a  +  (1  -  a.) Wn  (2.28) 

where  the  subscript  op  stands  for  “output”. 

A  high  initial  Kalman  gain  is  insured  by  setting  the  initial  conditions  as 


follows: 


for  i  =  1,2, . . .  N 


XIq|o[1,  i]  —  zo 


XI0|o[3,  i]  =  0 


Po|o[3, 3]  =  0 


z0 

(2.29) 

-  *o) 

(2.30) 

0 

(2.31) 

o\  x  10 

(2.32) 

T 2 

(2.33) 

Po|o[2, 1]  = 

(2.34) 

_  10<72 

Polo  [3, 1]  —  j,2 

(2.35) 

0 

(2.36) 

(2.37) 

\fU..rfL ..ttL  *T-4  A*.  it~- af V.~ . aft*. j/!«.  .* . «rL  at. 


3.  Estimation  Factors 


Some  factors  in  the  estimation  process  emerged  as  dominant  when  the 
observer  performance  was  evaluated  over  a  variety  of  conditions.  These  factors  are 
reviewed  in  this  section. 

a.  Hypothesized  Commands 

The  use  of  a  conditional  mean  between  the  set  of  hypothesized  com¬ 
mand  parameters  in  the  estimator  leads  to  the  dependence  of  the  estimate  on 
the  validity  and  accuracy  of  the  hypotheses  and  the  nature  (distribution)  of  the 
measurement  noise.  This  effect  can  be  demonstrated  by  the  following  simplified 
example.  Consider  an  unknown  scalar  parameter  U  which  is  to  be  estimated  from 
a  single  observation  z  =  U  +  n  where  n  is  a  random  variable  with  zero  mean  normal 
distribution  n  ~  (0, a2).  If  one  assumes  a  discrete  model  for  U  with  two  equally 
likely  possible  hypotheses  171  and  U 2  and  U2  >  f/l,  then  the  Bayes  mean-square 
estimate  is  given  by  the  conditional  mean 

u  =  E  {U\z}  p(ui I2)  =  L*7‘-  -  — V  P(U,)  (2-38) 

,  Jz 


where  /  is  the  probability  density  function  of  z.  This  estimate  is  shown  in 
appendix  B  to  yield 


U  =  UC  +  ^-tanh 


(2.39) 


where  Uc  and  U3  axe  the  center  and  separation  of  the  hypothesis  command  bank 
defined  by 


Uc 


Ul  +u2 

2 


U.  =  U2-  Ui 


(2.40) 

(2.41) 
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Y  OBSERVATION  [M/SECJ 
Fig.  2.3.  Conditional  mean  estimate. 

The  dependence  of  the  estimates  on  Uc  and  U,  and  on  the  noise  distribution  var: 
ance  a  is  clearly  seen  in  Eq.  (2.39).  The  function  U(z)  defined  by  Eq.  (2.39)  i 
plotted  in  Fig.  2.3.  The  bias  £  =  E  —  U  j  of  the  estimator  is  given  by: 

£  =  E  jz7c  +  y-tanh  J^|-(Z7  +  n  -  Uc)  -  U j  (2.42 

which  after  rearranging  and  extracting  the  constants  from  the  expectation  gives 

£  =Vc-U+^E{tmh£l(U  +  n-V')\.  (2.43 

Given  that  Uc  =  U  the  conditioned  bias  can  be  shown  to  be  zero  as  follows.  Sub 
stituting  U  —  Uc  into  Eq.  (2.43)  gives  the  bias  as 


£  =  ^r-E  ( tanh  ( ~t^\  \  . 


(0  44 


.-w  .  »  --  ■.  \  *-  v  a  v\  %  v  *.  v 

Since  the  tanh(x)  function  is  odd  and  the  density  function  of  n  is  even,  the  expec¬ 
tation  in  Eq.  (2.43)  is  zero. 

However  the  bias  is  non-zero  when  Uc  U,  i.e.,  when  the  hypothesis 
bank  is  not  centered  around  the  actual  value.  This  can  be  realized  for  example  by 
investigating  Fig.  2.3  for  values  of  U  larger  then  Ui  or  smaller  th*  .  U\ .  If  Uc  U 
then  U  will  be  incorrectly  estimated  with  the  amount  of  error  depending  on  the 
deviation  U  —  Uc. 

i 

This  simple  analysis  which  can  be  extended  to  the  MM  case,  provides 
|  a  new  insight  into  the  operation  of  the  MM  estimator  and  its  bias.  If  the  command 

\  bank  is  to  cover  all  possible  commands  of  a  target  which  could  move  at  full  speed 

i 

|  both  incoming  (decreasing  range)  and  outgoing  (increasing  range),  its  center  should 

i 

I  be  set  to  zero  speed.  While  a  choice  of  Uc  =  0  will  produce  good  results  when 

averaged  over  all  possible  target  scenarios,  it  will  produce  a  bias  E  for  every  specific 
scenario  with  average  command  E{U}  ^  Uc.  As  long  as  the  target  is  moving  in 
1  one  direction  (inward  or  outward)  the  center  of  the  hypothesis  bank  will  be  very 

[  different  from  the  average  command.  This  will  give  rise  to  an  estimation  bias. 

I 

In  the  example  shown  in  Fig.  2.4  the  command  U  of  a  target  driven  by 
a  true  command  of  10  m/sec  is  estimated  using  a  bank  of  7  filters  evenly  spanning 

l 

;  the  range  —15  to  +15  m/sec  centered  around  0  m/sec.  The  resulting  steady  state 

l 

|  command  bias  is  0.8  m/sec  as  shown  in  Fig.  2.4a.  The  position  error  resulting 

!  from  the  command  bias  is  150  m  (Fig.  2.4b). 
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b.  Basic  State  Separation 


The  dependence  of  the  MM  performance  on  the  hypotheses  separa¬ 
tion  warrants  some  further  investigation.  The  measurement  and  the  evaluation 
of  the  conditional  mean  Eire  performed  in  the  position  x(z)  domain.  The  system 
dynamics  as  represented  in  the  transition  matrix  and  the  Kalman  filter  translate  a 
hypothesized  command  to  a  hypothesized  position.  In  particular,  consider  a  system 


Xn  =^X„_,  +TU  +  *wn-l 


with  an  observation 


zn  —  HXn  4  Vfi 


modeled  by 


X„  =  *Xn_!  +  TUh  4- 


(2.45) 


(2.46) 


(2.47) 


Uh  represents  the  hypothesized  constant  command  here  assumed  to  be  different 
than  the  actual  constant  command  U.  The  observation  noise  v  and  process  noise  w 
are  both  normally  distributed  with  v  ~  N( 0,  cr£)  and  w  ~  N( 0,  a^,).  The  variables 
<f>,  r  and  H  are  the  transition,  control  gain,  and  observation  matrices  respectively. 
The  goal  is  to  find  the  mean  steady  state  deviation  £  of  X  from  X,  i.e., 

£  =  Jim^-E  {Xn|n  -  X„}  (2.48) 

resulting  from  the  mismatch  between  the  model  (U  h)  and  the  actual  system  ( U ). 
The  Kalman  prediction  is 

X„|„_i  =  ^Xn_i|„_i  +  TV  h  (2.49) 

and  its  average  is  given  by 


E  {*„!„_!  }  =  +E  {xn-l]n-i  }  +  TUh. 


(2.50) 


The  Kalman  estimate  Xnin  is 


X„in  =  (I  —  GnH)  Xnin_ i  -I-  Gnzr 


(2.51) 


which  after  substitution  of  Eq.  (2.45)  and  (2.46)  becomes 


Xn)n  =  (I  -  GnH)  ti,,.!  +  GnH  (^Xn_!  +  T U  +  Vwn-i )  +  Gnvn.  (2.52) 


Taking  the  expectation  of  Eq.  (2.52)  yields 

E  {xn|n  }  =  (I  -  G„H)  E  {x„|„-i }  +  GH*E  {Xn_j }  +  GnHH7 
+  GHtf£  K.J  +  G.,%} 
and  since  E{wn }  =  E{vn}  =  0  Eq.  (2.53)  becomes 


(2.53) 


E  jxn,n  }  =  (I  -  G„H)  E  {V* }  +  GnH^E  {Xn_x }  +  GnHTt/.  (2.54) 


If  Eq.  (2.50)  is  substituted  into  Eq.  (2.54)  the  result  is: 


E  {X„|„  }  =  (I  -  G„H)  (*E  }  +  r U„)  +  G„H(*£{X„-, }  +VU). 


(2.55) 


Now  taking  the  expectation  of  both  sides  of  Eq.  (2.45)  gives 


E  {X„}  =  $E  {X„_i}  4-  TU  *E  {u?n_j } 


(2.56) 


and  since  £{u>n}  =  0, 


£{Xw}ME{Xf,_1}+r£f. 


Subtracting  Eq.  (2.57)  from  Eq.  (2.55),  rearranging  terms  and  substituting  £„ 


defined  as 


£n  =  E  jxn|fl  -  Xn  j  =  E  |x„|„|  -  E  {Xn} 


(2.58) 
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results  in 


F  {£„}  =  (I  -  GnH)  *F  {£ n_i }  +  (I  -  G„H)  r (I7„  -  17) .  (2.59) 

At  the  steady  state  we  can  substitute  £  =  for  both  £n  and  £n-i 
as  well  as  G  =  G^  for  Gn,  in  Eq.  (2.59)  and  therefore  obtain 

£  =(I-GH)4£ +  (I-GH)T(Uh -U)  (2.60) 


or  alternatively 


[I  -  (I  -  GH)  4]  £  =  (I  -  GH)  T(UH  -  17).  (2.61) 

The  desired  transformation  of  the  speed  command  hypothesis  deviation  Lth  —  U 
to  the  steady  state  average  state  deviation  is  thus  found  to  be 

£  =  [I  -  (I  -  GH)^]"1  (I  -  GH)T(Uh  -  17)  (2.62) 

which  has  the  simple  form 

£  =  F (UH  -  U )  (2.63) 

where 

F  =  [I  -  (I  -  GH)  fl-1  (I  -  GH)  T.  (2.64) 

In  the  MM  the  hypothesized  commands  in  the  bank  deviate  from  the 
actual  command  value.  This  command  deviation  translates  into  the  deviation  of  the 
corresponding  states  from  the  actual  state  according  to  Eq.  (2.63).  The  separation 
of  the  states  in  the  bank  is  therefore  dependent  on  the  command  separation,  the 
system  matrices  and  the  steady  state  Kalman  gain. 


The  dominance  of  the  Kalman  gain  in  Eq.  (2.64)  is  clearly  evident.  A 
simple  scalar  example  will  demonstrate  the  effect.  Consider  an  attempt  to  estimate 
a  slowly  varying  constant  x  defined  by  the  model 

xn  —  ax  „_i  +  U  +  (2.65) 

where  w  is  a  process  noise  tv  ~  N  (0,ct^)  and  a  has  a  positive  value  close  to  but 
less  than  1.  The  observation  is  defined  by 

zn  =  Hxn  +  vn  (2.66) 

where  H  =  1  and  v  is  the  observation  noise  v  ~  N  (0,  er^).  If  the  model 

*n|n-l  =  ain- l|n-l  +Ujf+  Wn-1  (2-67) 

is  used  for  the  Kalman  filter  estimate  then,  the  model  mismatch  Uh  —U  will  give 
rise  to  a  steady  state  position  deviation.  Substituting  ^  =  a,r  =  l,H  =  l  into  Eq. 
(2.62)  gives  the  steady  state  position  deviation 

E  =  -  ~9  r  (U„  -  U)  (2.68) 

1  _(1  -g)a 

where  g  is  the  positive  scadax  steady  state  Kalman  gain  which  depends  on  the  vari¬ 
ances  ratio  of  the  process  and  the  measurement  noise.  The  deviation,  as  indicated 
by  Eq.  (2.68),  will  vary  from  a  minimum  of  0  (for  the  highest  gain  case  of  g  —  1) 
to  a  maximum  of  (  f°r  the  lowest  gain  case  of  g  =  0). 

For  the  MM  the  T,  H  matrices  and  the  allowable  range  of  maneuver 
commands  are  assumed  known.  Under  such  conditions  the  only  control  the  MM 
designer  has  over  the  state  separation  is  the  Kalman  gain.  A  high  gain  gives  heavy 
weight  to  the  new  measurement  which  is  common  to  all  the  filters;  this  tends  to 
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keep  the  filters  states  close  together  since  their  differences  are  not  emphasized.  A 
low  gain  gives  higher  weight  to  the  different  individual  channel  predictions  thus 
allowing  the  states  to  grow  apart.  The  states  diverge  until  their  average  innovation 
is  large  enough  to  compensate  via  the  Kalman  gain  for  the  command  mismatch. 

An  example  is  given  in  Fig.  2.5a  showing  the  positions  of  the  two 
extremes  and  the  center  filters  in  the  bank.  The  Kalman  gain  (position  component) 
was  G(l)  =  0.073  when  the  position  deviation  of  the  first  filter  was  1.4km.  Fig. 
2.5b  shows  a  similar  case  with  higher  gain  G(l)  =  0.041.  Note  that  the  separation 
between  the  first  and  center  tracks  is  increased  to  2.9Km  (note  the  different  scale). 


a.  G(l)=0.073 


b.G(l)  =  0.041 
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Fig.  2.5.  State  position  separation. 


In  Fig.  2.6  a  3-D  plot  showing  the  weight  vector  W  as  a  function  of 
time.  A  sharp  and  well  defined  peak  (a)  leads  to  a  less  noisy  and  overall  better 
maneuver  tracking  performance.  This  is  the  case  if  most  of  the  measurements  z 
(Eq.  (2.13))  fall  within  the  interval  spanned  by  the  state  bank  XI,  and  if  the 
innovation  variance  o\  is  of  the  order  of  the  position  separation.  If  many  of  the 
measurements  fall  outside  the  interval  spanned  by  XI,  as  is  the  case  when  the  gain 
is  too  high,  then  a  less  well  defined  peak  is  formed.  Such  a  case  is  shown  in  Fig. 
2.6b  where  the  resulting  “hesitation”  and  noisy  nature  of  the  weight  vector  may 
eventually  lead  to  loss  of  track. 

On  the  basis  of  the  forgoing  discussion,  it  can  be  seen  that  in  order  to 
minimize  the  estimation  bias  one  should  take  the  following  actions.  First  one  has 
to  ensure  that  the  state  bank  spans  the  expected  spread  of  position  measurements. 
This  can  be  done  by  tuning  the  Kalman  gain.  Secondly  one  should  maintain  the 
center  of  the  hypothesis  bank  close  to  the  actual  command.  This  can  be  done  by 
adaptively  relocating  the  command  bank  center  Uc  defined  as: 

1  N 

u'  =  if'Lu‘  <2-69> 

i=  1 

around  the  estimated  command  value.  Adaptive  recentering  has  not  been  at¬ 
tempted  in  previous  work  with  the  MM  tracking  filter. 

Trimming  the  Kalman  gain,  is  done  by  means  of  an  optimized  correc¬ 
tion  factor  Du.  Athans  and  Chang  [4]  refer  to  the  optimization  as  “more  of  an  art 
then  a  science.”  and  Ref.  9  attempts  to  derive  an  analytical  expression  for  Du*. 

*  The  argument  made  there  is  that  each  model  in  the  bank  is  subject  to  an 
additional  command  uncertainty  in  the  interval  around  the  discrete  value. 

This  is  of  limited  validity  since  the  command  uncertainty  for  every  channel  is  over 
the  entire  interval  Umaz  —  Um{n-  However  the  Singer  noise  component  which  is 
added  to  the  Du  compensates  by  maintaining  a  high  enough  Kalman  gain. 
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Fig.  2.6.  Weight  vector  history. 

Gain  control  by  means  of  the  Du  factor  reduces  the  need  for  the  Singer  state  com¬ 
ponent  w' ,  which  was  devised  partially  for  the  same  purpose.  Our  approach,  which 
was  found  to  be  more  effective  is  discussed  in  Section  D.  Adaptively  recentering 
the  command  bank,  was  a  more  involved  procedure  than  it  seem  to  be  at  first.  The 
method  is  also  described  in  Section  D. 


To  complete  the  discussion  Eq.  (2.63)  is  now  extended  to  case  of  a 
complete  MM  state  bank  with  N  channels.  Recall  from  Table  2.2  that  the  state 


bank  matrix  has  the  form 

XI  =  [X,|X2|...|Xtf]  (2.69) 

where  XiX2  ...X^  are  the  state  vectors  for  the  individual  filters,  and  that  the 
command  bank  vector  has  the  form 

UI  =  [UUU2,...UN]  (2.70) 

where  (7j,  U2,  . . .  Us  are  the  corresponding  commands.  The  state  bank  deviation 
matrix  defined  as 


£1  =  [£x ,  £2, . . .  £n]  =  [Xt  -  X,  X2  -  X, . . .  XN  -  X]  (2.71) 


is  given  by 


£  I  =  XI  -  X  N 


where 


N  =  [1,1,...  1). 


It  follows  from  Eq.  (2.63)  that  the  error  for  the  iih  filter  is  given  by 


£■>  =F  (U,  —  U) 


and  thus  £\  is  also  given  by 


£1  =  [F(Ut  -  U),  F(U2  -U),...,  F(Un  -  V )] 
£1  =  F(UI  -  U  ■  N) 


(2.72) 


(2.73) 


(2.74) 


(2.75) 

(2.75) 


namely 

XI  -  X  N  =  F(UI  -  UN).  (2.76) 

Now  if  each  control  U{  is  changed  by  the  same  amount  A U ,  then  the  change  in  the 
entire  command  bank  AUI  is  A U  •  N  and  the  new  command  vector  will  be  given 
by  UI  +  A U  •  N.  From  Eq.  (2.63)  the  corresponding  change  in  the  state  bank 
matrix  is  given  by 

AXI  =  F  ■  AUI  =  F  •  A£7  •  N.  (2.77) 

The  foregoing  discussion  implies  that  for  an  evenly  spaced  command  the  state 
bank  matrix  is  also  evenly  spaced  along  the  position  speed  and  acceleration  axis, 
the  state  spacing  AX  is  given  by 

AX  =  F  •  AU.  (2.78) 

Eq.  (2.78)  compared  very  favorably  with  simulation  results  (the  predicted  deviation 
was  within  ±0.01%  of  the  simulation  result.)  and  turned  out  to  be  very  useful  in 
ridding  the  MM  tracker  from  one  of  its  inherent  biases  (  in  Section  D). 

c.  Smoother  Order 

A  first  order  averager,  Eq.  (2.26),  (2.27),  and  (2.28)  was  added  to 
smooth  the  filter  output  in  order  to  reduce  the  output  variance.  The  price  paid  for 
the  smoothing  is  a  tracking  lag  error  for  targets  of  constant  speed  [9].  This  results 
from  the  mismatch  between  the  target  model  (second  order)  and  the  averager  (first 
order).  The  steady  state  emr  E  can  be  derived  using  the  final  value  theorem 

Xoc  =  lim  - - X(z)  (2.79) 

*->1  z 

where  X{z)  is  the  2  transform  of  xn.  Consider  the  scalar,  constant  rate  target 
position  xn  given  by 

xn  =  U0  ■  nu(n).  (2. SO) 


66 


where  u(n)  is  the  unit  step  function,  and  U0  is  the  rate  parameter  [position  units/ sample 
time  T],  and  a  smoothing  filter  of  the  form 


Vn  =  ayn-i  +  (1  -  a)xn.  (2.81) 

The  observation  z  is  here  assumed  to  be  noise  free  so  that  zn  =  xn •  The 
smoothing  filter  produces  a  lag  error  E  [position  units]  defined  as 

En  =  yn  —  xn.  (2.82) 


whose  value  at  steady  state  is  sought.  The  analysis  in  the  Z  domain  proceeds  as 
follows.  Let  'H(z)  be  the  system  response  and  X(z),y(z)  be  the  Z  transformations 
of  xn  and  yn.  Then 


•*(.>  =  -*  •  x  lu°  ■  u<-)]  =  U°z 


n.)  = 


dz 
1  —  az~x 
1  —  a 


{z  -  iy 


e(z)  =  y{z)  -  x{l)  =  (H{z)  -  1)X{Z)  = 


U0a  1 
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Applying  Eq.  (2.79)  yields  the  steady  state  error 
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and  the  steady  state  lag  error  is  thus 
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(2.83) 

(2.84) 

(2.85) 


(2.86) 


(2.87) 


indicating  large  errors  for  heavy  smoothing  (a  close  to  1). 

In  the  time  domain  Eq.  (2.81)  averages  the  past  estimated  value 
X (n  —  1)  with  the  current  measurement  of  x(n).  This  produces  an  error  even 
for  the  case  of  no  measurement  noise  since  the  position  x  is  different  at  the  two 
successive  time  samples.  This  intuitive  understanding  of  the  source  of  the  error  led 
to  the  development  of  the  improved  advancing  smoother  discussed  in  Section  D. 


d.  Coordinate  Linearization  and  Decoupling 


An  important  aspect  of  the  coordinate  decoupling  and  linearization  is 
that  the  target  helmsman  commands  in  the  decoupled  system  are  attached  to  a 
moving  coordinate  system,  namely  the  target  range  and  cross  range.  The  drawback 
of  such  a  system  is  that  a  (nonmaneuvering)  target  moving  with  constant  velocity 
along  any  straight  line  other  than  the  line  of  sight  has  to  be  described  as  maneu¬ 
vering  since  its  range  and  cross-range  speeds  are  constantly  changing.  However  in 
the  MM  the  probability  that  the  command  will  not  change  is  assumed  to  be  very 
high  (Pr  close  to  1). 

These  two  opposing  assumptions  produce  an  error  which  is  more  pro¬ 
nounced  at  short  ranges  where  range  and  cross  range  accelerations  are  higher.  At 
short  ranges  the  normalized  range  rate  (range  rate/range)  is  also  large.  This  im¬ 
plies  that  there  is  a  need  to  update  the  bearing  channel  transition  matrix  more 
frequently.  Under  the  foregoing  condition  the  advantages  of  this  linearization  and 
decoupling  approach  are  marginal. 

A  potential  remedy  to  the  coordinate  linearization  problem  may  be 
to  return  to  an  XY  Cartesian  system  or  to  adaptively  modify  the  probability 
transition  matrix  <j>  to  reflect  a  predicted  linear  motion.  These  ideas  were  not 
pursued  further  however  in  order  to  concentrate  effort  on  the  IH  effects  of  the 
MP  tracking. 

e.  Measurement  Noise 

In  passive  localization,  TDOA  from  the  target  via  different  acoustical 
paths  to  separate  parts  of  a  sonar  array  are  converted  to  depth,  range  and  bearing 
information.  The  process  is  error  prone  both  in  the  time  delay  estimation  and  in 
the  transformation  to  target  position  (range  depth  bearing).  This  is  one  of  the 
major  topics  addressed  in  this  research  and  most  of  Chapter  Three  is  devoted  to 


this  problem.  However  the  noise  dependence  on  range  affects  the  evaluation  of  the 
target  tracker  itself,  and  so  is  described  briefly  here. 

The  time  delay  estimation  error  depends  on  the  signal  to  noise  ratio 
(SNR)  of  the  received  signal.  The  SNR  in  turn  depends  on  the  transmission  loss 
and  therefore  indirectly  on  the  range.  Thus,  the  delay  noise  is  Tange  dependent* 
In  addition,  the  time  delay  estimation  limits  the  estimated  delay  to 
positive  values.  An  improved  noise  model  is  used  for  the  tracker  performance 
evaluation.  This  noise  model  is  described  in  the  Section  D. 

D.  IMPROVEMENTS  AND  MODIFICATION 

The  following  modifications  of  the  tracker  were  introduced  to  improve  it  and 
evaluate  its  performance  with  respect  to  the  issues  discussed  in  Section  C. 

•  An  adaptive  instead  of  a  fixed  hypothesis  bank  was  incorporated. 

•  A  second  instead  of  a  third  order  system  was  incorporated  in  the  model.  The 
second  order  system  was  demonstrated  as  sufficient  for  the  target  when  the 
steady  state  Kalman  gain  is  adjusted  by  the  correction  factor  Du  (see  eq. 
(2.16),  (2.17),  (2.20)). 

•  A  new  second  order  smoothing  algorithm  was  developed  and  used  to  replace 
the  first  order  smoother. 

The  modified  tracker  was  than  evaluated  with 

•  A  more  realistic  range  dependent  and  nonnegative  noise. 

•  Emphasis  on  range  bearing  coordinate  decoupling  and  linearization  at  short 
ranges. 


*  With  this  limitation  in  mind  the  disadvantage  of  the  coordinate  linearization, 
discussed  above,  at  short  ranges  is  further  emphasized. 
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a.  Concept 

The  problem  addressed  by  this  modification  is  the  bias  error  arising 
from  a  hypotheses  bank  which  is  not  symmetric  around  the  actual  command.  A 
possible  solution  to  this  problem,  cited  earlier,  is  to  adaptively  recenter  the  hypoth¬ 
esis  bank  around  the  actual  value.  Recall  that  when  the  center  Uc  of  the  command 
bank  (referred  to  here  as  “the  center”  for  convenience)  is  equal  to  the  value  of 
the  true  command,  this  bias  is  removed.  Since  the  actual  value  of  the  command 
is  not  known,  our  algorithm  uses  the  command  estimate  Uop  instead.  The  main 
idea  is  therefore  to  periodically  recenter  the  command  bank  around  the  estimated 
command  value. 

A  straightforward  implementation  is  to  average  the  estimated  control 
over  a  time  interval  long  compared  to  the  system  and  the  observer  time  constants 
and  to  feed  back  the  average  as  the  new  center.  Specifically  the  difference  between 
the  estimated  command  U0p  and  the  center  Uc  is  averaged  using  an  autoregressive 
filter  to  produce  the  center  deviation  Ucd- 


U cd  n  —  ^U cd  n—  1  (1  o)(£^op  n  —  1  Uc  n) 


(2.88) 


This  average  value  is  used  to  dynamically  shift  the  center  and  the  complete  com¬ 
mand  bank. 


Uc  n-fl  —  U c  n  "1"  ^  cd  n 

UIn+1  =  UI„  ■+■  Ucd  n 


(2.89) 

(2.90) 


Although  this  method  was  implemented,  its  performance  was  poor  for  the  following 
reason.  The  combination  of  the  command  hypothesis  UI  and  the  weight  vector 
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W  is  a  discrete  description  of  the  command  probability  distribution  estimated  by 
the  Bayesian  recursion  (Eq.  (2.23)).  Updating  the  command  hypothesis  vector  but 
keeping  the  weight  vector  unchanged  effectively  modifies  the  command  distribution 
and  disrupts  the  recursive  command  estimation. 

Further,  the  command  hank  yields,  over  some  time,  a  corresponding 
bank  of  hypothesized  states  in  the  form  of  the  matrix  XI  which  based  on  Eq.  (2.76) 
is 

XI  =  \£XE2  . .  •  £n|  +  X  •  N  =  F(UI  +  U  •  N). 


Different  states  correspond  to  the  new  updated  command  bark.  Shifting  the  com¬ 
mand  hypotheses  bank  introduces  an  undesirable  transient  in  the  state  vectors  of 
the  filters.  This  transient  seriously  degrades  the  tracking.  Furthermore,  the  feed¬ 
back  of  the  state  vector’s  transient  into  the  recursive  command  estimation  process 
(Eq.  (2.22)  and  (2.23a)),  may  lead  to  instability  of  the  entire  tracker. 

It  is  clear  therefore  that  updating  the  command  bank  center  alone 
disrupts  the  recursive  estimation  process  unless  there  is  a  corresponding  update 
of  the  weight  vector  W  and  the  state  matrix  XI.  In  simulation  recentering  the 
command  bank  without  updating  the  weight  vector  and  the  state  matrix  resulted 
in  complete  divergence  of  the  command  estimate, 
b.  Complete  Model  Updating 

An  alternative  approach  was  devised  where  the  complete  model  was 
updated  and  the  instability  was  eliminated.  The  method  was  as  follows: 

•  The  update  was  restricted  to  occur  at  discrete  events  in  time  when  the  averaged 
deviation  of  the  center  Uc  from  the  estimate  Uop  reached  a  preset  threshold. 

•  The  commands  and  the  command  center  were  restricted  to  take  on  a  set  of 
discrete  quantized  equally-spaced  values.  Updating  was  done  by  shifting  the 
entire  bank  up  or  down  one  position. 
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•  The  complete  model  including  command  and  state  banks,  and  the  weight 
vector  were  updated. 


The  time  discretization  was  intended  to  prevent  the  continuous  dy¬ 
namic  feedback  which  led  to  instability.  The  time  constant  used  for  the  exponential 
command  deviation  averaging  was  set  larger  than  both  the  system  and  the  MM 
estimator  time  constants.  This  was  done  to  ensure  complete  transient  recovery 
from  one  update  before  another  update  is  allowed. 

The  model  update  was  implemented  in  a  cyclic  manner.  The  least 
likely  filter  in  the  bank,  i.e.,  the  one  corresponding  to  the  command  furthest  away 
from  the  estimated  command,  is  dropped.  A  new  and  more  likely  command  is 


added  at  the  closer  end  of  the  bank.  The  probability  of  the  new  channel  is  initially 

N 

set  to  that  of  the  one  that  is  dropped,  so  that  the  total  probability  ]T]  Wt  is  kept 

t=i 


equal  to  1.  If  the  update  takes  place  at  time  k  then 


UIfc+  =  UIfc_  ±  At/  •  N 


(2.91) 


where  the  subscripts  k  and  represent  the  control  prior  to  and  after  the  update 
and  the  sign  of  At/  depends  on  the  sign  of  the  averaged  center  deviation  Ucd  (Eq. 
(2.88)). 

A  block  diagram  of  the  modified  MM  is  shown  in  Fig.  2.7  and  Fig. 
2.8  shows  an  example  of  a  model  update  process  assuming  the  update  takes  place 
at  time  index  n  =  k.  The  variables  at  time  k~ ,  i.e.,  just  prior  to  the  update, 
are  shown  in  Fig.  2.8a;  the  variables  immediately  after  the  update  (time  k+ )  are 
shown  in  Fig.  2.8b.  Prior  to  the  update  the  bank  is  centered  around  Uc  —  0 
auad  spans  the  interval  —15  to  +15  [m/sec]  as  shown  along  the  horizontal  axis. 
The  second  (lower)  horizontal  axis  represents  the  first  component  (position)  of 
the  state  vectors  of  all  filters  in  the  bank  (this  is  the  first  row  of  XI)  before  the 


update.  This  spans  the  interval  8.5  to  11.5  km.  During  the  update  the  first 
hypothesis  (f/i  )k-  =—15  m/sec  is  dropped  and  a  new  hypothesis  {U-r  )k+  is  added 
at  +20  m/sec  with  a  corresponding  state  position  XI[1,  7]*+  at  12  km  (Fig.  2.8b). 
The  weights  (conditional  probabilities)  for  all  the  filters  except  that  of  the  new 
one  are  shifted  one  position.  This  is  done  in  order  to  maintain  the  correspondence 
with  the  hypotheses.  The  new  filter  assumes  the  weight  of  the  filter  that  is  dropped 
(0.01  in  the  example  shown).  The  command  bank  is  now  properly  centered  around 
5  m/sec.  The  command  estimation  transient  is  shown  in  Appendix  C  to  be  given 
by 

U transient  =  NAI7(W1)*_  (2.92) 

For  the  specific  case  in  Fig.  2.8  this  transient  is  therefore 

U transient  =  7  •  5  •  0.01  =  0.35m/sec.  . 

This  minor  transient  is  further  smoothed  by  the  output  smoother  and  its  effect  on 
the  output  U0?  is  negligible. 

The  command  in  the  bank  and  the  center  Uc  are  restricted  to  the 
evenly  spaced  discrete  values.  The  center  update  is  further  limited  to  a  one  position 
shift  up  or  down  from  its  current  position  on  the  discretized  scale.  The  command 
bank  update  is  thus  implemented  by  means  of  the  linear  shift  given  by  Eq.  (2.91). 
Since  the  states  of  the  bank  are  also  evenly  spread  over  the  state  scales  (position, 
speed  and  acceleration).  The  corresponding  update  of  the  stafes  matrix  XI  can 
thus  be  implemented  by  means  of  a  linear  shift  as  follows 

XI*+  =  XIt_  ±  AX  •  N  (2.93) 

where  AX  is  the  constant  vector  F  ■  AU  ■  N  (Eq.  (2.78)).  The  complete  model  is 
thus  updated  by  a  simple  one  position  shift.  The  shift  is  linear  for  the  command 
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Fig.  2.S.  Complete  model  update. 
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i  UI  and  the  state  XI,  and  circular  for  the  weight  vector  W.  The  simplicity  of 

i  the  model  update  procedure  (which  requires  no  computation)  resulted  from  the 

constraint  of  the  center  Uc  to  assume  values  on  the  discrete  command  scale. 

Tracking  with  the  complete  adaptive  model  (updating  UI,  W  and  XI) 
is  shown  in  Fig.  2.9.  Notice  the  command  and  position  biases  developing  just  prior 
to  the  update  2.9c  which  is  reduced  by  the  two  updates  at  the  18<A  and  35th 
minute,  and  the  almost  unnoticeable  transient.  Further  persistence  of  the  same 
target  motion  would  lead  to  additional  updates  which  would  eventually  reduce  the 
bias  completely. 

2.  Second  Order  Target  Model 

In  Eq.  (2.8)  the  second  order  target  model  was  augmented  to  a  third  order 
model  in  order  to  account  for  the  uncertainty  in  the  maneuver  command.  Prac¬ 
tically  this  augmentation  effected  the  tracking  by  increasing  the  error  covariance 
matrix  P  which  in  turn  increases  the  Kalman  gain  (Eq.  (2.14)).  Since  the  Kalman 
gain  is  so  critical  in  the  MM  estimator  it  is  additionally  trimmed  by  the  factor  Du. 

In  the  original  third  order  Singer  model  the  third  state  variable  w  was 
the  only  state  variable  that  accounted  for  the  target  acceleration.  In  the  MM 
the  command  (and  the  corresponding  acceleration)  is  separately  estimated  as  the 
independent  state  U.  With  both  the  acceleration  estimate  and  the  Kalman  gain 
handled  otherwise  the  use  of  the  third  state  w  is  redundant  and  can  be  eliminated. 
Elimination  of  this  state  results  in  the  following  second  order  model 
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XI2  n  =  n  — 1  +  r2Un_i 


(2.94) 


Where  XI2H2  and  T2  are  the  first  two  rows  of  XI,  H  and  T  respectively,  and 
is  the  2x2  upper  left  block  part  of  ^  (see  Eq.  (2.15)).  This  configuration  will  save 
computation  time  and  simplify  the  gain  adjustment  procedure. 
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Fig.  2.9.  Tracking  of  an  updating  model 


Fig.  2.10  compares  the  tracking  of  a  third  order  model  with  steady  state 
Kalman  position  gain  G[l]  =  0.041  to  the  tracking  of  a  second  order  model  with 
practically  the  same  position  gain  (G[l]  =  0.042).  The  gain  of  the  second  order 
model  is  controlled  by  <xu  which  affects  the  gain  via  D„  =  IV„r7\  The  performance 
of  the  two  models  is  similar  except  the  third  order  system  has  a  larger  overshoot. 

3.  Second  Order  Smoother 

An  improved  second  order  smoother  referred  to  as  an  advancing  smoother 
was  introduced  [25]  in  order  to  reduce  the  constant  speed  target  lag  error  (Eq. 
(2.75)).  The  basic  idea  is  to  average  the  new  filter  output  xn|n  with  the  predicted 
point  xn|n_j  rather  than  with  the  previous  estimated  position  xn_!|n_i.  The 
current  estimated  speed  is  used  to  generate  the  predicted  position  according  to 

(Xop)n  =  a  •  [1,  T,  0]  •  (Xop)„_,  +  (1  -  a)[l,  0, 0]X[l]n.  (2.95) 

If  carried  out  to  the  full  extent,  this  principle  would  turn  the  smoother  into  a 
Kalman  filter  by  itself.  The  need  for  this  additional  smoothing,  which  is  usually 
handled  by  means  of  a  reduced  Kalman  gain,  results  from  the  optimizing  the  gains 
to  meet  the  overall  MM  requirements.  This  forces  the  additional  external  smooth¬ 
ing. 

4.  Improved  Measurement  Noise  Model 

A  range  dependent,  nonnegative  noise  is  used  in  this  simulation  in  place 
of  the  range  independent  normally  distributed  noise  used  in  past  work. 

a.  Dependence  on  SNR 

The  dependence  of  the  delay  estimation  noise  standard  deviation  (STD ), 
on  the  range,  via  its  dependence  on  the  SNR.  was  modeled  here  in  two  ways.  In 
the  first  model  the  range  dependence  is  lumped  into  the  following  equation. 

STD,(R)  =  STDlo  Rp  (2.96) 
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where  STDia  is  the  i-th  TDOA  error  initial  standard  deviation  at  500  m  from  the 
source,  and  p  is  the  delay  noise  range  power  factor  which  depends  on  propagation 
loss.  The  variable  p  takes  on  values  in  the  range  of  1  to  3.  A  low  value  for  p  (p  =  1 ) 
corresponds  to  a  high  SNR  and  cylindrical  propagation  while  a  higher  value  (p  =  2) 
corresponds  to  low  SNR  and  spherical  propagation.  Higher  p  values  corresponding 
to  increased  loss  due  to  ray  bending  were  not  considered. 

In  the  second  model,  a  more  recent  and  detailed  description  was  used 
based  on  Ref.  14.  Here  the  exact  dependence  of  the  variance  on  the  bandwidth, 
the  observation  time  T  and  the  SNR  was  implemented.  The  parameters  of  the  first 
simplified  noise  model  (ST Do  and  p)  were  set  to  provide  an  overall  similar  noise 
variance  which  is  similar  to  the  results  of  the  second  noise  model, 
b.  Nonnegative  Time  Delay 

Another  feature  of  the  time  delay  estimation  for  a  single  receiver  is 
the  inability  to  distinguish  a  leading  signal  replica  from  a  lagging  one.  This  results 
from  the  symmetry  of  the  autocorrelation  function  (ACF)  of  real  signals.  The 
nonnegative  time  delay  effect  was  simulated  by  applying  the  absolute  value  operator 
to  the  time  delay  after  the  noise  was  added.  That  is 


n  m  =  k,  + 1\|  (2.97) 

where  v,  ~  JV(0,  STD^(R)).  This  operation  produces  a  bias  primarily  for 
stem  R)  <  1’  which  is  also  discussed  in  Appendix  D. 

After  the  addition  of  the  noise  the  simulated  TDOA  undergoes  a  non¬ 
linear  inversion  to  depth  and  range  which  also  produces  a  bias.  The  bias  was 
studied  for  the  idealized,  homogeneous  and  space  invariant  delay  noise,  by  Moose 
in  Ref.  11.  Investigation  of  this  error  and  its  correction  for  the  IH  case,  and  the 
means  to  account  for  it  are  discussed  in  Chapter  Four. 
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c.  Depth  Range  and  Bearing  Noise 

When  the  filter  was  evaluated  by  itself  (without  the  prefilter)  noise 
was  added  to  the  depth  and  range  measurement  directly.  Normally  distributed 
range  dependent  noise  was  used  such  that 


Dm  =  D  +  vd 
Rm  =  R  +  vr 


(2.90) 

(2.99) 


where 


vi  ~N(0,STDl(R)) 
vr  ~  N(0,  STDi(R)) 


(2.100(a)) 

(2.101(6)) 


and  the  standard  deviations  given  by 


STDi(R)  =  STDdo  ■  Rp> 
STDr(R)  =  STDro  ■  RVi 


(2.101(a)) 

(2.101(6)) 


with  ST Di0  STDro  being  the  initial  STD  at  500  m  and  pb  the  range  dependence 


factor. 


Range  dependent  noise  was  added  to  the  bearing  channel  as  well  in 


order  to  maintain  a  consistent  uniform  simulation.  The  bearing  measurement  is 


where 


Bm  =  B  +  vb 


Vb  ~  N(0,  STDb). 


(2.102) 


(2.103) 
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The  bearing  error  standard  deviation  used  here  is  given  by: 

STD(R)b  =  STDbo  ■  Rpb  (2.104) 

where  STDf,0  is  the  bearing  error  standard  deviation  at  500  m  from  the  source. 

E.  SIMULATION 

1.  Model  Description 

The  simulation  model  was  designed  to  evaluate  the  overall  MP  tracker 
performance,  both  in  the  idealized  homogeneous  case  and  later  in  the  realistic  and 
IH  environment  (Chapter  Four).  The  following  were  set  as  subgoals  for  the  design: 

•  Support  of  a  multilevel  system  evaluation  starting  from  the  individual  algo¬ 
rithm  up  to  the  complete  integrated  system. 

•  Provide  a  detailed  evaluation  of  the  3-D  maneuvering  target  tracker  perfor¬ 
mance  with  emphasis  on  the  3  axis  integration  and  the  tracker  modifications. 

•  Support  of  a  gradual  departure  from  the  idealized  homogeneous  MP  config¬ 
uration  towards  the  more  realistic  IH  case.  The  transition  includes  both  the 
simulated  IH  measurement  and  the  algorithm  designed  to  deal  with  this  dis¬ 
tortion.  This  again  is  in  support  of  Chapter  Three. 

In  order  to  maintain  a  consistent  and  unified  simulation  environment  the 
model  was  designed  to  meet  all  of  the  above  subgoals  by  optional  use  of  its  various 
components.  The  model  is  discussed  here  in  its  entirety;  however  the  details  of 
some  of  its  components  are  addressed  in  the  following  chapters. 

A  block  diagram  of  the  simulation  model  is  shown  in  Fig.  2.11.  The 
simulation  proceeds  as  follows.  Depth,  range  and  bearing  target  position  D ,  R ,  B 
is  generated  by  the  scenario  generator.  The  depth  and  range  are  converted  to 
TDOA  rj ,  r2  by  the  MEDIUM.  Noise  generated  by  the  NOISE  module  is  added  to 
the  TDOA  to  form  the  measured  rlmr2rn  which  are  fed  to  the  PREFILTER.  The 
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prefilter  performs  a  functional  inversion  to  compute  the  depth  and  range  measure¬ 
ments  Dm  Rm.  Noise  is  also  added  to  the  bearing  to  simulate  a  realistic  bearing 
measurement  Bm.  The  TARGET  TRACKER  uses  the  three  noisy  measurements 
to  produce  a  3-D  target  estimate.  The  estimate  includes  position  Rop,Bop,  rate 
RoP,  Bop  and  command  Ur  op,  U&  op,  f or  the  range  and  cross  range  axis,  and  posi¬ 
tion  Dop  for  the  depth  axis.  The  noisy  measurements  and  the  tracker  estimates 
are  compared  to  the  output  of  the  scenario  generator  to  produce  the  measurement 
and  estimation  errors.  Details  of  each  of  the  blocks  are  described  below. 

The  SCENARIO  GENERATOR  contains  a  model  for  both  the  target  and 
its  pilot.  It  includes  the  target  dynamics  (Eq.  (2.10))  as  well  as  a  prescribed  target 
trajectory  and  it  outputs  the  target  position  in  depth,  range,  and  bearing.  Three 
different  types  of  pilots  were  programmed. 

Pilot  1  describes  the  target  motion  by  an  initial  state  and  a  set  of  discrete 
and  independent  sequences  of  maneuvering  commands  along  each  of  the 
three  axis. 

Pilot  2  is  similar  to  Pilot  1  except  it  provides  a  continuous  gradual  change  of 
the  maneuvering  command. 

Pilot  3  prescribes  a  target  moving  along  a  straight  horizontal  line  with  constant 
speed  (Eq.  (2.8)).  As  discussed  in  Section  C  this  target  will  be  seen  as 
maneuvering  by  a  cylindrical  coordinate  tracker. 

The  MEDIUM  converts  the  target  depth  and  range  coordinates  to  time 
delay  differences  T\ ,  r2  (or  ti ,  t2).  It  includes  a  model  for  the  effects  of  the  medium, 
the  receiver,  and  the  noise.  Two  groups  of  models  were  programmed,  namely  the 
homogeneous  group  and  the  inhomogeneous  group  described  below. 


The  homogeneous  group. 

Medium  1  This  is  a  bypass  mode  where  the  depth  and  range  are  used  directly, 
without  the  conversion  to  time  delays  and  back  to  depth  and  range. 
This  mode  is  used  to  evaluate  the  MM  tracker  by  itself. 

Medium  2  The  approximate  formulas  valid  for  R»  D,  [Ref.  9] 
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T2 


2D0(Dq  -  D) 

RC 

2  (Dw  —  Dq^Dv  —  Do  +  D) 

RC 


were  used  for  the  direct  function. 


(2.105a) 

(2.1056) 


Medium  3  The  exact  closed  form  analytical  expression  of  time  delays  as  func¬ 
tion  of  depth  and  range  derived  by  Hassab  (Eq.  (1.1))  were  used  here. 
The  implied  assumptions  here,  besides  straight  line  propagation,  are 
that  the  delays  are  perfectly  resolvable  and  can  be  associated  with  their 
corresponding  paths. 
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The  inhomogeneous  group. 

Medium  4  This  is  the  only  medium  that  was  implemented  in  the  group.  Here 
the  three  assumptions  implied  above  were  removed  and  an  IH  medium 
with  finite  resolution  and  non-associable  time  delays  is  generated.  The 
theory  and  implementation  of  Medium-4  are  the  subject  of  Chapter 
Three.  TDOAs  for  the  realistic  IH  medium  are  noted  as 


The  NOISE  models  were  used  for  the  TDOA  measurements.  They  were 
noise  type  (NT)  one  and  two.  NT  1  used  Eqs.  (2.96)  and  (2.97).  NT  2  used  Eq. 
(2.96)  and  the  more  elaborate  noise  model  which  provides  the  STD,  as  a  function 
of  range  and  the  target  signal  to  noise  ratio  at  range  of  1  m  (SNR0).  The  model 
and  the  parameters  it  uses  are  described  in  Appendix  D. 

When  no  medium  and  prefilters  were  used,  noise  was  added  directly  to 
depth  and  range  using  Eq.  (2.98)  and  (2.99).  Noise  was  also  added  to  the  bearing 
according  to  Eq.  (2.102). 
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The  PREFILTER  converts  the  noisy  TDOA  estimates  transforming  them 
to  depth  and  range.  This  is  done  by  inverting  the  nonlinear  direct  function.  Six 
different  prefilters  were  designed  and  grouped  into  two  groups:  homogeneous  and 


inhomogeneous  prefilters. 

Homogeneous  prefilters. 

Prefilter  -  1  This  is  the  bypass  case  where  no  functional  inversion  is  done  since 
the  inputs  are  depth  and  range,  corresponding  to  medium  1. 


Prefilter  -  2  This  uses  the  following  approximated  relation 

D  =  Co - - 

r‘  + 

2DqDw 


C  (T>  +  d^d;t 0 


(2.106a) 


(2.1066) 


valid  for  R  »  D  corresponding  to  medium-2  [Ref.  9], 

Prefilter  -  3  This  is  the  exact  analytic  closed  form  relation  of  Eq.  (1.2)  corre¬ 
sponding  to  Medium  3. 

Inhomogeneous  prefilters. 

Prefilters  -  4,5,6  These  correspond  to  the  inhomogeneous  medium,  with  finite 
TDOA  resolution  and  non  associable  delays  of  medium  4.  These  again 
are  discussed  further  in  Chapter  Three. 


The  TRACKER  uses  the  measured  depth,  range  and  bearing  as  inputs  for 
the  3-D  maneuvering  MM  target  tracker.  Third  and  second  order  target  models 
with  first  and  second  order  smoothers  were  realized.  The  command  bank  was  either 
fixed  or  adaptive  as  discussed  in  Section  D.l. 

2.  Description  of  the  Output  Plots 

Three  main  plots  aire  produced  ais  output  for  each  axis  of  the  cylindrical 
coordinate  system;  they  are  the  position,  the  speed  and  the  command  plots.  When 
a  combination  of  plots  is  presented  the  plot  symbol  markings  are  synchronized  in 


time  on  all  the  plots  of  a  given  run  to  provide  a  means  for  cross  reference  between 
them. 

The  values  of  the  simulated  Rims  were  recorded  and  plotted  only  every 
fifth  sample  time  to  avoid  cluttering  of  the  picture.  An  example  of  a  full  sample 
rate  plot  (with  a  point  plotted  for  every  sample  time  period)  versus  a  reduced 
sample  rate  (point  plotted  every  fifth  sample  time)  is  shown  on  Fig.  2.12. 

The  position  plots  (depth,  range,  bearing,  and  X,  Y  as  functions  of  time) 
include  the  true,  measured,  and  estimated  positions.  The  true  position  is  the 
scenario  generator  output  D,R,B ;  the  measured  position  is  the  prefilter  depth 
and  range  outputs  Dm ,  Rm  and  noisy  bearing  Bm ;  and  the  estimated  position  is 
the  tracker  output  Dop  Rop  Bop.  Position  error  plots  show  the  difference  between 
the  measured  or  estimated  positions  and  the  true  position. 

The  speed  plots  present  range  or  cross  range  speed  as  a  function  of  time. 
Various  combinations  of  the  actual  and  estimated  speed  and  speed  command  are 
included.  The  center  of  the  command  bank  is  also  shown.  Note  that  when  Pilot-3 
was  used,  constant  Ux  and  U9  are  applied  (Eq.  (2.8)).  The  commands  Ur  and  Ub 
axe  not  used.  The  model  however  does  estimate  Ur  and  Ub,  which  Eire  presented 
on  the  plots. 

Command  weights  of  the  recursive  Bayesian  estimator  axe  plotted  in  3-D 
plots  of  the  weight  vector  eis  a  function  of  time.  These  plots  provide  a  good  measure 
of  the  MM  Eidjustments  and  tracking  quality. 

3.  Parameters  Selection 

The  large  number  of  parameters  required  for  each  run  dictated  a  strict  and 
well  controlled  mechanism  of  parameter  selection.  This  was  especially  required  for 
the  realistic  IH  CEise  where  there  was  a  large  amount  of  ocean  and  acoustic  data 
and  parEimeters.  The  data  was  divided  into  the  following  three  categories:  target 
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and  scenario;  medium  and  prefilter;  MM  parameters.  The  detailed  definition  of 
the  parameters  and  their  values  in  the  simulation  is  presented  in  Appendix  D. 

F.  SIMULATION  RESULTS 

Five  straight  line  trajectories  are  presented  differing  mainly  in  the  noise  used. 

Run  1  used  range  dependent  depth  and  range  noise  without  the  multipath 
measurement  ( Mdn  =  1,  Pfn  =  1).  The  range  dependence  p  was  set  at  1  and  the 
depth  and  range  error  STD  at  500  m  were  set  to  1  and  50  respectively. 

Runs  2-5  used  range  dependent  TDOA  noise  based  on  the  improved  noise 
model  {NT  =  2)  with  SNR  and  p  set  according  to  Table  2.3. 

A  detailed  list  of  the  other  parameter  settings  is  included  in  Appendix  E. 

TABLE  2.3 

NOISE  PARAMETERS 

Run  SNRo  [d£]  p 

2  50  1 

3  70  2 

4  60  2 

5  50  2 

Results  of  Run  1  are  shown  in  Fig.  2.13.  Range  (a)  depth  (b)  and  range  error 
(c)  are  shown.  The  increase  of  measurement  and  tracking  errors  with  increasing 
range  is  seen  on  all  three  plots.  The  very  effective  noise  filtration  of  the  filter  is 
demonstrated  in  Fig.  2.13c  after  the  completion  of  the  initial  transient  (around 
the  12th  min  and  on).  The  range  measurement  noise  STD  which  was  1000  m  at  10 
km  was  reduced  to  less  than  100  m.  Range  error  reaching  17km  develops  around 
the  CPA  which  occurs  at  the  17,A  minute  at  a  range  of  3  km.  The  results  of  this 


RANGE  [KM] 


a.  RANGEACTUAL  •.  MEASURED  +  .  ESTIMATED  7 


b.  DEPTH  ACTUAL  •.  MEASURED  +  ,  ESTIMATED  V 


C.  RANGE  ERROR-.MEASURED  +  ESTIMATED  7 


TIME  TWIN] 

Fig.  2.13.  Run-1:  Range  dependent  DR  noise. 


Run  set  a  reference  performance  for  the  MM  maneuvering  target  tracker  by  itself 
(without  the  prefilter). 


Results  of  Rim  2  are  shown  in  Fig.  2.14.  The  performance  of  the  MM 
tracker  and  the  homogeneous  medium  and  prefilter  are  demonstrated  with  range 
dependent  TDOA  noise. 

As  can  be  seen,  the  cylindrical  propagation  (p  =  1)  provides  very  favorable 
conditions  even  when  the  initial  SNRo  is  relatively  low  (50  dB  at  1  m  range).  Note 
the  four  model  update  cycles  the  MM  goes  through  in  response  to  the  changing 
range  rate.  Also  note  the  vanishing  range  tracking  error  which  results  from  the 
use  of  the  advancing  smoother  ( 25tk  minute  in  Fig.  2.14a).  Good  bearing  tracking 
and  the  combined  horizontal  tracking  are  shown  for  this  run  in  Fig.  2.15.  The  XY 
tracking  error  at  the  CPA  is  clearly  evident. 

The  effect  of  increasing  TDOA  noise  as  a  result  of  decreasing  target  signal 
strength  (Run-3,  4)  is  shown  in  Figures  2.16  and  2.17  The  failure  of  the  tracker 
in  Run  4  past  the  range  of  15Km  at  the  27th  minute  of  this  run  is  also  evident  in 
the  command  estimation  plot  shown  in  Fig.  2.17c. 

Very  noisy  TDOA  measurements  resulting  from  low  SNRq  and  spherical 
spreading  (Run-5)  is  shown  on  Fig.  2.18.  Large  estimation  bias  which  results  from 
the  nonnegative  effect  is  clearly  seen  at  ranges  larger  than  12  km.  The  depth  is 
also  underestimated  at  these  ranges. 
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G.  SUMMARY  OF  TRACKER  EVALUATION 

The  performance  of  the  modified  target  tracker  is  summarized  below. 

•  Generally  the  tracker  is  well  suited  for  3-D  maneuvering  target  tracking. 

•  Measurement  noise  variance  is  significantly  reduced  by  the  tracker.  The  range 
measurement  STD  of  1000  m  can  be  easily  handled  and  reduced  to  around  50 


•  The  performance  depends  strongly  on  range,  especially  when  range  dependent 
MP  measurements  are  used.  However  up  to  15-17  km  the  filter  tracks  well 
without  readjustments. 

•  Tracking  errors  that  can  reach  2  km  in  range  develop  at  short  ranges  for 
nonmaneuvering  targets  movie  r  along  straight  lines  due  to  the  coordinate 
decoupling  and  linearization. 

•  Recentering  of  the  command  bank  around  the  average  actual  command  is 
required  in  order  to  remove  the  inherent  MM  estimation  bias.  The  new  model 
update  technique  devised  here  effectively  accomplishes  this  task  by  recentering 
the  command  bank  around  the  estimated  command  value. 
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AND  RANGE  MEASUREMENTS 


A.  INTRODUCTION 

The  following  idealizing  assumptions  mentioned  in  Chapter  One  and  used  in 
previous  broadband  MP  tracking  studies  are  removed  in  this  chapter: 

•  Homogeneous  straight  line  propagation. 

•  No  propagation  and  reflection  loss. 

•  Infinitely  resolvable  TDOA. 

•  Ability  to  associate  delays  with  acoustical  path. 

The  impact  of  these  effects  on  the  tracker  is  investigated  in  Section  B.  Section 
C  presents  the  improved  prefilter  which  provides  the  functional  inversion  of  time  to 
depth  and  range  measurement  for  the  IH  realistic  medium.  Chapter  Four  analyzes 
the  performance  of  the  new  inversion  algorithm  as  a  component  in  the  overall  MP 
tracking  system. 

B.  REALISTIC  MEDIUM  AND  RECEIVER  EFFECTS 

I.  Propagation  Loss  and  Reflections 

The  most  dominant  effect  the  medium  has  on  the  propagating  acoustic 
wave  is  the  transmission  loss.  Three  contributors  to  the  loss  axe  discussed:  the 
spreading,  the  absorption,  and  the  imperfect  reflections  from  the  boundaries.  Spe¬ 
cial  attention  is  given  to  the  effect  that  the  reflection  from  the  non-ideal  boundaries 
has  on  the  multipaths. 
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a.  Spreading  Loss 

The  spreading  loss,  as  the  name  implies,  is  the  result  of  the  spreading 
of  the  acoustic  energy  over  a  growing  area,  as  the  distance  from  the  source  increases. 
In  spherical  spreading,  the  energy  spreads  evenly  over  the  surface  of  a  sphere.  This 
makes  the  intensity,  the  power  per  unit  area,  inversely  proportional  to  the  squared 
range.  The  relation  is 

I(R)  =  /„  -  R~2  (3.1) 

where  I0  is  the  intensity  at  a  distance  of  1  m  from  the  source. 

Under  other  propagation  conditions  different  range  dependency  re¬ 
sults.  For  example,  with  cylindrical  propagation  the  intensity  is  inversely  propor¬ 
tional  to  the  range  itself,  that  is 

I(R)  =  Jo  •  fT1  (3.2) 

In  deep  water,  away  from  the  surface  or  the  bottom,  one  finds  spherical 
propagation  to  be  a  good  assumption  for  short  ranges.  In  shallow  water,  where 
MP  effects  are  more  dominant,  cylindrical  propagation  is  a  better  representation. 

In  the  IH  case  the  wavefront  is  distorted  by  the  variation  in  the  speed 
of  sound.  Although  the  resulting  spreading  loss  is  not  simple,  a  dependence  on 
the  p-th  power  of  range  can  still  be  useful  as  a  first  order  approximation.  The 
approximation  is 

I(R)  =  Jo  •  R'p  (3.3) 

where  the  empirical  constant  p,  which  we  call  the  range  propagation  power,  is  not 
limited  to  integers. 


b.  Absorption  Loss 

Part  of  the  energy  of  the  propagating  acoustic  wave  is  absorbed  by 
the  water  and  eventually  is  turned  into  heat.  The  absorption  loss  depends  on  the 
salinity  of  the  water,  the  frequency,  and  the  distance  traveled.  In  short  to  medium 
ranges  (below  15  km),  and  at  sufficiently  low  frequency  (below  10  KHz)  this  is  not 
a  dominant  phenomenon  and  will  therefore  not  be  considered  here. 

c.  Reflections  From  the  Boundaries 

The  laws  of  pressure  and  velocity  continuity  across  the  boundary  gov¬ 
ern  the  reflections  from  the  surface  and  the  bottom  [Ref.  23].  A  wave  propagating 
in  medium  1  impinging  on  a  boundary  at  an  angle  0i  is  both  reflected  back  at  an 
angle  /?i  and  refracted  into  medium  2  at  an  angle  /?2  as  shown  in  Fig.  3.1. 


Fig.  3.1.  Reflection  from  a  boundary. 


The  acoustic  pressure  reflection  coefficient  defined  as 


r=  — 


where  pi,pr  are  the  acoustic  pressures  of  the  incident  and  reflected  waves  in  medium 
1.  The  reflection  coefficient  is  given  by 

r  =  (3.5) 

Zi  cos  (3\  +  Z\  cos  02 

where  Z\  and  Z2  are  the  specific  acoustic  impedances  in  media  1  and  2,  respec¬ 
tively.  The  specific  acoustic  impedance  is  defined  as  the  ratio  between  the  complex 
amplitude  of  the  acoustic  pressure  p  and  that  of  the  magnitude  of  the  particle 
velocity  U.  The  angle  $2  is  given  by  $2  =  sin-1  sin  where  C\  and  C2  are 
the  speeds  of  sound  in  the  two  media. 

The  specific  acoustic  impedance  for  a  plane  wave  can  be  expressed  in 
terms  of  the  speed  of  sound  C  and  the  equilibrium  density  p0  of  the  medium  as 


Z=po-C 


which  is  known  as  the  characteristic  impedance. 

Since  the  characteristic  impedance  of  air  is  much  smaller  than  that  of 
water,  the  surface  boundary  has  a  reflection  coefficient  close  to  —1  (see  Eq.  (3.5)). 
Further,  since  the  characteristic  impedance  of  the  bottom  is  generally  larger  than 
that  of  the  water,  the  resultant  reflection  coefficient  of  the  bottom  is  positive. 

The  reflections  from  the  water  air  boundary  becomes  somewhat  more 
involved  when  the  layer  of  air  bubbles  formed  near  the  surface  by  the  wind  and 
the  ocean  waves  is  considered.  But  unlike  the  air  (which  is  isotropic)  the  bottom 
may  be  lossy  and  anisotropic.  In  particular,  shear  waves  excited  by  the  impinging 
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acoustic  wave  travel  in  the  bottom  with  different  characteristics  than  the  refracted 
longitudinal  waves.  These  effects  give  rise  to  less  predictable  and  complex  acous¬ 
tic  impedance  and  complex  reflection  coefficients  which  strongly  depend  on  the 
bottom  type.  For  a  boundary  with  frequency  dependent  characteristic  impedance 
Z2  =  r2  +  jx 2,  the  pressure  reflection  coefficient  will  vary  both  with  frequency  and 
angle  of  incidence,  and  will  depend  on  the  bottom  type  and  sea  state  (for  surface 
reflection).  Exact  prediction  of  the  phase  of  a  reflected  wave  is  made  difficult  by  the 
large  variation  in  pressure  reflection  coefficients.  This  phenomena  has  a  marked 
influence  on  MP  tracking.  The  intensity,  which  is  related  to  the  square  of  the 
pressure,  has  the  reflection  coefficient  T/  given  by 

r  /  =  |if  (3.7) 

which  also  depends  on  the  above  ocean  and  acoustic  wave  parameters.  In  a  very 
simple  practical  analysis,  the  energy  loss  due  to  the  various  boundary  effects  is 
lumped  into  a  single  number  called  the  reflection  loss.  The  reflection  loss  depends 
mostly  on  the  sea  state,  the  bottom  type,  and  the  frequency  band.  Loss  values  in 
the  range  of  3  to  30  dB  per  bounce  are  commonly  found.  The  sensitivity  of  the 
reflection  phase  to  frequency  and  angle  of  incidence  has  a  marked  influence  on  MP 
measurements. 

2.  Inhomogeneous  Medium 

a.  Speed  of  Sound  Variations 

The  speed  of  sound  in  the  water  is  given  by  the  formula  [23] 

C  =  1449.2  +  4.6T  -  0.055 T  +  0.00029T  +  (1.34  -  0.01  T)(S  -  35)  +  0.016Z?  (3.8) 


where  T  is  temperature  in  degrees  centigrade  5  is  the  salinity  in  parts  per  thousand 
and  D  is  the  depth  in  meters. 
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The  sound  speed  is  obviously  not  constant  in  the  vertical  water  column, 
due  to  change  of  the  pressure  and  the  temperature  with  depth.  The  variation  of 
sound  velocity  with  depth,  referred  to  as  the  sound  velocity  profile  (SVP)  is  the 
main  source  of  inhomogeneity  in  the  acoustic  ocean  medium.  Waves  propagating 
in  such  a  medium  are  refracted  in  a  complicated  manner,  and  obey  the  second 
order  linear  partial  differential  scalar  wave  equation  [20] 

=  Xm(f,r)  (3.9) 

where  <p(t,  r )  is  the  velocity  potential  at  time  t  and  position  r.  C(r)  is  the  speed  of 
sound  at  position  r,  and  xm(t,  r)  is  the  source  distribution  (volume  flow  rate  per 
unit  volume).  The  acoustic  particle  velocity  vector  U(<,  r)  and  pressure  p(t,  r )  axe 


given  by 


pO.r)  =  -po^O.r) 
U(t,r)  =  V<?(t,r). 


(3.10) 

(3.11) 


Solution  of  the  wave  equation  for  the  general  case  is  very  difficult. 
However,  when  the  speed  of  sound  is  a  linear  function  of  depth  the  approximate 
ray  acoustic  solution  yields  simple  closed-form  equations.  The  detailed  develop¬ 
ment  of  the  ray  acoustics  approximation  is  based  on  the  WKB  approximation  and 
is  presented,  for  example,  in  Ref.  20.  The  resulting  ray  tracing  equations  are 
presented  below. 

b.  Ray  Tracing 

The  parameters  and  geometry  of  a  typical  ray-trace  are  shown  in  Fig. 
3.2.  Note  that  the  rate  of  change  of  sound  speed  with  depth  is  a  constant  known 
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Fig.  3.2.  Acoustic  ray  path. 

as  the  sound  speed  gradient  g ,  which  has  units  of  [sec-1],  and  the  speed  of  sound 


C  is  given  by 


C(D)  =  C0+gD 


(3.12) 


where  C0  is  the  speed  of  sound  at  the  surface. 

In  transition  from  depth  Di  where  the  sound  speed  is  C\  to  depth  D2 
with  sound  speed  of  C2 ,  the  ray  angles  0\  and  02  obey  Snell’s  law 


sin  0\  _  Ci 
sin  fc  C2 


(3.13) 


Thus  a  ray  emanating  from  a  source  at  D\,R\,  at  an  angle  0\ ,  reaches  depth  D2 


at  an  angle  02  given  by 


02  =  sin" 


■\  ( C"!  '  a\  •  -1  ^0 

(cT  =  s'„  [cf 


+  2  .  „ 

— n  smft 
+  gD  1 


(3.14) 


and  at  a  range  R2  given  by 


R2  =  Ri  - ~~r  (cos 0i  -  cos02) 


(3.15) 


which  can  be  evaluated  by  solving  Eq.  (3.14)  for  02  and  substituting  into  eq. 
(3.15).  The  travel  time  T  and  the  arc  length  S  of  the  path  are  given  by 


T  =  1]ntan(^/2) 
g  tan(/?i/2) 


(3.16) 


S--fV(A-A). 


(3.17) 


The  ray  traverses  the  depth  range  (DR)  plane  along  of  a  circle  with  radius  Rr 


given  by 


r>  Co  Cl  COS  0i 

Kray  —  — 


(3.18) 


centered  at  the  initial  range  and  the  conceptual  depth  Dctnttr-  This  is  the  point 
where  the  speed  of  sound  would  have  become  zero  if  the  gradient  had  been  constant 
to  that  point;  it’s  value  measured  from  the  surface  is  therefore  given  by 


'center  — 


(3.19) 


Tracing  the  ray  intensity  is  based  on  conservation  of  energy.  A  hypo¬ 
thetical  amount  of  radiated  power  trapped  in  a  ray  tube  produces  intensity  which 
is  inversely  proportional  to  the  cross  sectioned  area  of  the  tube.  Tracing  the  cross 
sectional  area  along  the  tube,  using  the  ray  tracing  equations,  provides  the  required 


intensity  computation. 


Aoi.fli)  _  S(Di,R7) 
7(d,,r,)  S \Dur,) 


(3.20) 


A  singular  case  arises  when  the  cross  sectional  area  of  the  ray  tube 
reduces  to  zero  and  drives  the  intensity  to  infinity.  This  is  called  a  focal  point. 
Focal  points  in  which  the  intensity  is  very  high  are  actually  encountered  in  the  ocean 
and  the  acoustic  intensities  measured  there  are  very  high.  Ray  acoustic  intensity 
computations  are  obviously  not  valid  under  such  conditions.  A  ray  going  through 


a  focal  point  is  phase  shifted  by  90  deg.  This  phenomenon  which  is  experimentally 
measurable,  can  be  explained  analytically  although  not  in  a  trivial  manner. 

Ray  equations  allow  tracing  a  ray  in  a  medium  of  constant  SVP  gra¬ 
dient,  however,  a  constant  SVP  gradient  is  rarely  the  case  in  practice.  When  the 
SVP  is  not  linear,  it  can  be  represented  by  a  piecewise  linear  function  which  has 
constant  gradients  within  depth  layers.  Eq.  (3.14)  through  (3.17)  are  then  used  to 
trace  a  ray  inside  a  given  layer.  The  continuity  of  the  linearized  SVP  and  Snell’s 
law  ensure  ray  angle  continuity  in  the  transition  between  layers.  Consecutive  trac¬ 
ing  of  the  ray  inside  the  constant  gradient  layers  can  thus  construct  the  complete 
ray  path  throughout  the  entire  medium, 
c.  Effect  on  MP 

(1)  Time  Delay.  An  example  will  help  demonstrate  the  effect  of 
the  IH  propagation  on  MP  depth  and  range  measurement.  A  set  of  time  delays 
for  am  IH  case  is  computed  for  a  given  target  and  set  of  oceam  conditions  using 
Eq.  (3.16).  These  time  delays  are  then  used  as  independent  variables  of  the 
homogeneous  inverse  function  Eq.  (1.2).  The  erroneous  inverted  depth  and  range 
are  compared  to  the  original  values. 

A  constant  gradient  of  g=  0.05  sec-1  in  the  100  m  neair  the 
surface, *in  water  depth  of  513  m  and  a  target  at  depth  of  100  m  and  range  of 
7  km,**  produce  TDOAs  to  a  receiver  at  depth  of  162  m  of  tx  =  5.48  ms  and 
r2  =  28.92  ms.  When  substituted  into  Eq.  (1.2),  which  is  based  on  straight  line 
propagation,  a  depth  of  19  m  and  range  of  5.64  km  are  produced.  This  represents 
a  percent  error  of  20%  in  range  relative  to  the  true  7  km  range. 


(2)  Lack  of  Closed  Form  Solution  -  The  Eigenray  Problem. 

Of  significant  importance  is  the  fact  that  in  all  the  tracing  equations  the  initial 
transmission  angle  at  the  source  (ft)  is  assumed  to  be  known.  This  angle,  which 
can  alternatively  be  replaced  by  the  reception  angle  (ft),  actually  ‘selects’  the  ray 
to  be  traced. 

The  reverse  problem  of  finding  the  terminal  angles  of  a  ray  pass¬ 
ing  through  two  given  end  points  is  called  the  eigenray  problem.  The  eigenray 
problem  does  not  have  a  closed  form  solution  in  the  IH  case,  especially  not  in  the 
multilayered  SVP  case.  This  is  part  of  the  difficulty  which  prevented  compensation 
for  the  medium  inhomogeneity  in  previous  work  and  thus  became  a  key  challenge 
of  our  work. 

(3)  Loss  of  Direct  Ray.  Another  outcome  of  the  IH  propagation 
is  the  possible  loss  of  direct  path  between  the  source  and  receiver  due  to  the  ray 
bending.  Fig.  1.3a,  which  was  computed  by  a  Fortran  program  [26]  using  the  ray 
tracing  equations,  shows  an  example  of  a  receiver  Rx  placed  in  an  area  not  reached 
by  direct  rays  from  source  Tx  which  is  referred  to  as  a  “shadow  zone.”  The  first 
arrival  in  this  case  is  the  reflection  from  the  bottom.  The  second  arrival  is  the  ray 
reflected  from  the  surface  and  then  from  the  bottom.  The  third  arrival  is  the  one 
reflected  from  the  bottom  first  and  then  from  the  surface.  The  above  mentioned 
difficulty  to  predict  the  exact  amplitude  and  phase  of  the  reflected  wave  is  further 
enhanced  for  a  case  of  double  reflection  such  as  the  one  shown.  The  polarity  of 
the  ACF  peaks  becomes  less  predictable.  Thus,  it  is  obvious,  that  the  IH  medium 
amplifies  the  problem  of  path  identification. 

In  Comparing  Fig.  1.3a  to  Fig.  1.1,  the  difference  in  the  MP 
structure  is  mainly  that  the  direct  path  T0  of  Fig.  1.1  does  not  exist  in  Fig.  1.3b.  A 
scalar  called  the  bounce  count  (BC)  was  devised  here  in  an  attempt  to  characterize 
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the  MP  structure.  It  is  defined  as  the  total  number  of  depth  extremes  in  the  set 
of  the  first  three  resolvable  arrivals  (two  TDOAs). 

The  bounce  count  of  the  straight  line  propagation  of  Fig.  1.1  is 
2  due  to  the  single  bounces  of  both  the  surface  and  the  bottom  paths.  The  bounce 
count  of  the  first  three  arrivals  of  the  MP  structure  in  Fig.  1.3b  is  4  due  to  the 
single  bounces  of  the  first  and  second  arrivals  and  the  two  bounces  of  the  third 
path.  The  bounce  count  therefore  distinguishes  between  the  two  MP  structures. 

While  ambiguous  as  an  absolute  descriptor,  the  bounce  count 
was  found  to  be  very  useful  in  practice  primarily  in  identifying  a  change  in  the 
MP  structure.  By  plotting  the  BC  over  a  depth  range  cross  section  of  the  ocean, 
it  became  easy  to  divide  the  plane  into  regions  of  consistent  propagation  as  will 
be  shown  later.  Care  should  be  exercised,  however,  in  detailed  analysis  using  the 
bounce  count  since  radically  different  MP  structures  could  produce  identical  bounce 
counts. 

3.  Ocean  Data 

The  travel  time  of  paths  reflected  from  the  bottom  is  obviously  dependent 
on  the  ocean  depth.  Error  in  the  assumed  depth  will  cause  errors  in  the  inverted 
target  depth  and  range.  The  CRLB  for  the  joint  estimation  of  ocean  depth  and 
target  depth  and  range  was  developed  by  Hamilton  [Ref.  19]  as  a  composite  estima¬ 
tion  problem,  with  depth  assumed  normally  distributed.  The  impact  of  a  constant 
depth  error  is  demonstrated  in  Chapter  Four. 


The  effect  of  the  limited  delay  resolution,  and  the  inability  to  associate 
delays  with  paths  are  considered  here  along  with  an  outline  of  a  hypothetical  but 
realistic  form  of  delay  estimation  instrumentation. 


a.  Bandwidth  and  Delay  Resolution 


! 

1 
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The  target  broadband  noise  has  a  limited  bandwidth.  Further,  the 
high  frequency  portion  of  the  spectrum  tends  to  attenuate  in  a  shorter  distance 
than  the  lower  frequency  portion  as  the  sound  travels  through  the  water,  due  to 
absorption  processes.  The  receiver,  on  the  other  hand  is  constrained  by  engineering 
design  compromises  arising  primarily  from  the  beamformer,  and  computational 
load  associated  with  its  sampling  rates.  The  result  is  that  the  acoustic  receiver 
bandwidth  is  rarely  ever  more  then  a  few  KHz  wide. 

The  limited  bandwidth  translates  into  a  limited  resolution  between 
adjacent  time  lags.  If  one  assumes  a  uniform  signal  spectral  density  over  the  entire 
frequency  band  of  0  to  5  Hz,  the  main  lobe  of  the  autocorrelation  function  of  the 
original  signal  z(t)  will  be  of  the  form 

ACFzz{r)  =  ACFZZ( 0)  •  sinc(r  •  5)  (3.21) 

and  its  first  zero  occurs  at  1/5  sec.  Replicas  of  the  original  signal  which  are  less 
then  1/5  apart  will  form  a  nonresolvable  combined  ACF  peak.  (In  the  extremely 
unlikely  event  of  two  opposite  sign  equal  amplitude  replicas  the  ACF  peak  may 
even  disappear  altogether.) 

b.  Nonidentiflable  Path 

The  combination  of  the  effects  discussed  earlier  render  the  association 
of  the  TDOA  with  the  corresponding  path  a  difficult  task  indeed  which  is  imprac¬ 
tical  in  many  target  tracking  situations.  The  polarity  of  the  ACF  peaks  is  strongly 
affected  by  the  following: 


I 


3 
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•  Potential  loss  of  the  direct  reference  path. 

•  Complex  frequency  and  angle  dependent  reflection  coefficients. 

•  Phase  shift  through  a  focal  point. 

An  alternative  approach  taken  here  is  to  rely  only  on  the  time  difference  values, 
and  to  disregard  the  ACF  polarity  and  its  possible  association  with  a  specific  path. 

The  response  of  the  time  delay  tracker  to  this  nonresolvable  nonasso- 
ciable  situation  will  determine  the  behavior  of  the  target  tracker.  For  this  reason, 
the  specific  instrumentation  of  the  time  delay  estimator  assumed  here  is  briefly 
described  below. 

5.  Instrumentation  Outline 

A  discrete  set  of  ACF  values  spanning  the  range  of  expected  lag  time,  is 
computed  by  time  or  frequency  domain  methods,  using  the  time  sampled  acoustic 
signal.  A  polynomial  is  fitted  to  the  sampled  ACF  lag  data  and  its  first  two 
amplitude  maxima,  are  located  (excluding  t  =  0).  The  lags  times  corresponding 
to  these  peaks  are  the  MP  time  differences  of  arrival  (or  delays).  The  delays  jure 
usually  further  filtered  by  a  simple  exponential  averager.* 

The  first  two  lags  are  preferred  for  the  inversion  since  they  result  from 
simpler  paths  and  are  expected  to  have  higher  SNR.  But  when  any  two  of  the 
first  three  ACF  peaks  coincide  only  one  combined  and  somewhat  erroneous  peak  is 
formed.  An  example  where  this  problem  is  significant  and  very  typical  is  a  source 
placed  near  the  surface.  The  direct  and  surface  reflected  paths  from  the  source 
to  a  distant  receiver  will  have  very  similar  path  length  and  travel  time,  and  their 
ACF  peak  will  not  be  resolvable. 

•  Another  technique  of  similar  consequence  uses  interpolation  between 
adjacent  ACF  samples  to  reconstruct  the  continuous  delay,  Friedlander  [16]. 
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Other  lags  corresponding  to  more  complex  and  longer  paths,  do  exist  and 
can  be  used  in  such  cases.  Having  longer  lag  time,  these  delays  will  not  vanish,  due 
to  loss  of  resolution,  incidentally  with  the  first  two.  The  third  or  higher  lags  can 
thus  provide  an  alternative  measurement,  at  least  while  one  of  the  first  two  lags  is 
lost. 

In  the  instrumentation  model  simulated  here  a  third  ACF  lag  is  tracked  all 
along,  and  is  substituted  as  the  second  lag  only  when  any  of  the  first  three  arrivals 
coincide.  The  two  nonresolvable  lags  are  reported  as  one  lag  under  this  situation. 
Coincidence  is  defined  with  some  safety  margin  that  ensures  substitution  of  the 
third  lag  prior  to  the  actual  loss  of  resolution. 

In  Fig.  3.3  plot  of  the  first  three  TDOA  as  a  function  of  history  time  k 
is  shown.  At  time  T#,  the  ACF  peak  of  rx  is  not  resolvable  from  the  main  lobe 
at  r  =  0,  and  is  consequently  replaced  by  rj  which  in  turn  is  replaced  by  r3.  The 
third  TDOA  r3  was  tracked  but  not  output  before  ?V2.  If  rj  and  r2  merge  at  THX , 
as  shown  in  Fig.  3.4  their  combined  peak  will  be  reported  as  rx  and  r3  will  replace 
T2- 

In  the  straight  line  homogeneous  propagation,  identification  of  the  mea¬ 
sured  delays  as  tx  and  r2  is  assumed  such  that  they  can  be  used  in  Eq.  (1.2).  In  the 
realistic  IH  case  the  delays  are  identified  instead  as  the  first  (shorter)  and  second 
(longer)  lags.  A  vector  of  two  time  delays,  the  first  of  which  is  always  the  smaller, 
is  output  by  the  estimator  and  these  modified  TDOAs  are  denoted  as  tx  and  r2 . 

6.  Resulting  Ambiguity 

The  proposed  method  of  time  delay  extraction  immediately  raises  the  ques¬ 


tion  of  ambiguity.  There  can  obviously  be  more  then  one  position  with  the  same 
two  time  delays  as  demonstrated  by  the  simple  following  case. 
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For  an  observer  at  depth  162  m  and  water  depth  503  m  the  targets: 


Position 

Target  1 

Target  2 

Depth  m 

-100 

-242 

Range  m 

6000 

4782 

produce  To ,  T\ ,  T2  direct  surface  and  bottom 

travel  time,  respectively  of 

Travel  time  [sec] 

Target  1 

Target  2 

To 

4.00055 

3.192079 

T\ 

4.00997 

3.210297 

t2 

4.01877 

3.201499 

yielding  travel  time 

differences  of 

TDOA  [msec] 

Target  1 

Target  2 

II 

1 

9.41 

18.2 

3 

II 

1 

£ 

18.2 

9.41 

After  sorting  these 

become 

Modified  TDOA  [msec]  Target  1 

Target  2 

U  =  Minin,  t2 } 

9.41 

9.91 

#2  =  Max{rx,  r2} 

18.2 

18.2 

Yielding  the  same  n  t 2  set  for  two  vastly  different  depth  range  positions. 


This  ambiguity  is  a  direct  result  of  the  inability  to  associate  paths  and 
delays.  Means  to  reduce  this  ambiguity  to  a  manageable  proportion  are  covered  in 
Section  C.  A  more  effective  use  of  the  high  bounce  count  paths  designed  to  further 


reduce  the  ambiguity  is  discussed  in  Chapter  Six. 


Effects 

a.  Introduction 

The  2-D  direct  function  transformation  of  source  depth  and  range 
(DR)  position,  to  MP  time  delays,  (ri ,  T2),  can  be  viewed  as  surfaces  of  delays  over 
the  2-D  depth-range  position  plane.  These  surfaces  are  the  basis  of  the  MP  inver¬ 
sion  procedure  and  understanding  their  nature  is  important  in  the  development  of 
the  inversion  process. 

The  graphical  interpretation  of  the  direct  function  will  be  presented 
in  this  section  including  the  following  effects  which  the  medium  and  receiver  have 
upon  this  function: 

•  Lack  of  path  identification. 

•  Limited  resolution. 

•  IH  propagation. 

Four  types  of  surfaces  are  presented  starting  with  a  theoretical  case  of 
homogeneous  propagation,  and  resolvable  and  identifiable  multipaths.  The  above 
listed  effects  are  then  added  in  three  accumulating  steps  ending  with  the  realistic 
case  of  IH  propagation,  limited  resolution  and  non  identifiable  paths. 

b.  Case  Definition 

The  number  of  parameters  determining  the  MP  structure  is  very  large 
and  attempting  to  investigate  all  of  their  various  combinations  is  an  impossible 
and  not  a  valuable  effort.  A  number  of  representative  cases  were  selected  and 
the  parameters  involved  were  coded  into  cases.  A  case  name  is  composed  of  one 
letter  and  four  digits,  e.g.,  C1123,  which  has  the  following  meaning.  The  leading 
character  defines  the  type  of  assumptions  and  instrumentation  used.  The  first 
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digit  defines  the  SVP  and  the  second  digit  defines  the  ocean  and  observer  depth. 
The  last  two  digits  define  a  grid  of  depth  range  points  for  which  the  TDOAs  are 
computed.  A  detailed  description  and  settings  of  the  various  parameters  is  given 
in  Appendix  E. 

To  demonstrate  the  nature  of  the  direct  function  surfaces  three  cases 
were  used.  The  first  case  for  an  ocean  depth  around  500  m  and  receiver  depth  of  162 
m  (Ullll)  represents  homogeneous  propagation  and  infintely  resolvable  TDOAs 
which  are  associable  with  the  paths  for  bottom  depth  of  503  m  and  receiver  depth 
of  162  m.  The  second  case  (Si 111)  represents  the  same  case  but  without  the 
associability  of  TDOA  with  the  multipath.  The  third  case  (Cl 111)  represents 
the  same  case  with  both  the  TDOA-path  associability  and  the  perfect  resolution 
assumptions  removed.  Three  dimensional  (3-D)  surface  and  contour  plots  of  TDOA 
Ti ,  7*2  or  ,  t-2  over  the  DR  plane  are  plotted  for  each  case.  Regions  of  consistent 
MP  structure  are  indicated  by  contours  of  bounce  count  over  the  DR  plane.  A 
short  description  points  to  the  issues  of  interest  in  the  plots, 
c.  TDOA  Surfaces 

Homogeneous,  perfectly  resolvable  and  identifiable  Ullll. 
TDOA  surfaces  and  contour  lines  for  this  ideal  case  are  shown  in  Fig.  3.5. 

Both  surfaces  are  monotonic  and  T\  can  be  either  smaller  or  larger  then  r2.  Both 
delays  decrease  with  range  but  opposite  changes  in  rj  r2  result  when  depth  is  chang¬ 
ing  along  a  fixed  range.  Note  that  t\  is  small  and  r2  is  large  close  to  the  surface, 
and  the  reverse  is  true  near  the  bottom.  No  bounce  count  plots  are  presented  since 
the  count  is  2  all  over  the  plane.  Note  the  vanishing  t\  and  r2  near  the  surface 
and  the  bottom  respectively  which  causes  the  merging  of  the  corresponding  ACF 
peaks  with  the  main  ACF  lobe  at  r  =  0.  These  lags  are  substituted  by  those  from 
multiple  reflection  paths  in  the  limited  resolution  case.  Also  note  that  the  TDOA 


8000 


12000 


4000 


8000 


12000 


RANCE  [U] 

Fig.  3.5.  Ideal  TDOA  surface  (Ullll) 
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vanishes  with  range,  this  tends  to  make  the  effect  of  delay  estimation  noise  more 
significant  at  longer  ranges.  Since  in  addition  the  TDOA  estimation  noise  increases 
with  range,  the  MP  method  becomes  inherently  limited  to  short  ranges. 

Homogeneous  sorted  perfect  resolution  Sllll  modified  TDOA 
T\ ,  r2  surfaces  and  contours  for  this  ideal  case  where  association  of  the  path  and 
delays  is  not  assumed  are  shown  in  Fig.  3.6.  When  the  TDOAs  are  sorted  to 
produce  t\  and  <2  as  described  earlier  the  result  is  t\  <  t2  by  definition.  The 
monotonocity  of  the  surfaces  is  lost  along  the  DR  loci  where  t\  =  t2.  The  trend  of 
reduced  delays  T\ ,  r2  at  long  ranges  has  not  changed.  The  depth  still  determines 
the  relation  between  T\  and  r2  but  not  in  the  simple  unambiguous  manner  as  in 
the  unsorted  case. 

It  is  only  T\  now  that  vanishes  near  the  surface,  when  it  becomes 
smaller  than  the  delay  estimation  resolution  it  will  be  replaced  by  what  is  here 
plotted  as  r2.  The  r2  will  then  be  replaced  by  the  arrival  lag  of  the  next  multipath, 
the  one  bouncing  from  both  the  surface  and  the  bottom.  Since  this  will  only  happen 
for  the  realistic  receiver,  the  bounce  count  here  (perfect  resolution)  is  still  2  all  over 
the  DR  plane  and  is  therefore  not  plotted.  Contours  of  TDOA  in  the  first  100  m 
near  the  surface  are  shown  in  Fig.  3.7  as  a  reference  for  the  later  comparison  to  an 
IH  case. 

Homogeneous  realistic  finite  resolution  Cllll.  Fig.  3.8  shows 
the  resulting  surfaces  when  the  TDOA  are  not  associable  with  the  multipath 
(sorted)  and  a  finite  delay  resolution  (0.5  ms)  is  assumed. 

The  main  difference  between  this  and  Sllll  case  is  the  substitution 
of  a  high  order  reflection  time  delay  when  time  delays  from  simpler  reflections  are 
not  resolvable.  This  comes  about  near  the  surface  where  rj  is  small  and  at  depth 
of  around  -176  m  relative  to  the  receiver  where  ti  and  r2  are  similar  causing  the 
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replacement  of  T2  by  the  lag  corresponding  to  the  higher  bounce  count  path  (t3 ). 
The  bounce  count  plot  is  helpful  in  recognizing  this  effect  which  of  course  is  also 
obvious  on  the  TDOA  surface  and  contours  as  well. 


C.  IMPROVED  FILTER 
1.  Concept  Overview 

The  subject  of  this  section  is  inversion  of  realistic  IH  multipath  time  differ¬ 
ence  of  arrivals  (TDOA)  to  source  depth  and  range  coordinates.  Since  there  is  no 
closed  form  solution  for  either  the  direct  or  the  inverse  relation  between  time  delays 
and  position,  two-step  numerical  solution  is  developed.  As  before  the  function 
that  computes  depth  and  range  from  time  delays  is  referred  to  simply  as  the  direct 
function,  and  the  procedure  that  computes  depth  and  range  from  time  delays  is 
referred  to  as  the  inverse  function.  The  latter  may  at  times  be  multivalued  and  so 
may  not  be  a  one-to-one  function  in  the  mathematical  sense. 

First,  the  direct  function  values  of  the  TDOA  are  computed  over  a  dis¬ 
crete  grid  of  target  depth  and  range  positions.  This  computation  is  done  off-line. 
Measured  time  delays  are  then  inverted  to  obtain  position  (depth  and  range).  This 
computation  must  be  done  on-line. 

In  the  off  line  computation  an  eigenray  model  is  used  to  find  the  travel 
times  of  the  multiple  paths  between  the  receiver  and  a  given  depth  and  range 
point.  The  travel  time  produced  is  further  processed  to  simulate  the  realistic  delay 
estimation.  Generation  and  processing  of  the  TDOA  is  repeated  for  every  point  in 
the  entire  DR  grid.  The  processed  TDOA  data  thus  produced,  is  stored  as  a  table 
of  time  delay  as  a  function  of  the  depth  and  range. 
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In  the  on  line  step,  a  set  of  two  measured  time  delays  is  inverted  to  produce 
depth  (and  range)  using  2-D  linear  interpolation.  Depth  (and  raDge)  values  of 
three  points  in  the  direct  function  table,  which  have  time  delay  values  closest  to 
the  measured  set,  are  used  as  a  basis  for  the  interpolation.  To  locate  these  points 
a  unique  search  algorithm  is  developed.  The  algorithm  takes  advantage  of  a  priori 
information  available  in  the  form  of  the  predicted  target  position.  The  prediction, 
in  turn,  is  derived  from  past  depth  and  range  measurements  processed  by  the  target 
tracker. 

The  use  of  the  eigenray  model  and  the  delay  estimation  simulation  to 
generate  the  direct  function  is  described  in  Section  C.2.  The  inversion  process  is 
described  in  Section  C.4. 

2.  Direct  Function 

The  precomputed  tabulated  direct  function  of  the  modified  TDOA  over 
the  DR  grid  is  computed  by  means  of  an  eigenray  model  which  is  repetitively 
exercised  over  the  DR  grid  points.  This  process  is  described  next. 

a.  The  Eigenray  Model 

An  eigenray  model  is  a  program  that  finds  the  multiple  rays  which  path 
between  two  given  end  points,  providing  their  angles,  path  length  and  travel  time. 
Computing  the  TDOA  of  the  MP  structure  requires  very  accurate  computation 
of  the  path  travel  time.  This  results  from  the  TDOA  being  a  difference  of  a  few 
milliseconds  between  travel  times  which  are  of  the  order  of  a  few  seconds. 

In  the  search  for  an  eigenray  model  only  compact  programs  requiring 
small  computation  time  on  readily  available  microcomputers  were  considered.  This 
was  done  in  order  to  maintain  compatibility  both  with  the  computational  resources 
available  to  this  research  and  those  available  in  onboard  systems.  Since  no  closed 
form  solution  for  the  eigenray  problem  is  available  (recall  discussion  in  Section 


III.B.2.c)  the  eigenray  models  axe  iterative  in  nature.  The  conflicting  requirements 
for  high  accuracy,  limited  word  length  and  short  run  times  axe  met  by  only  a  few 
existing  models. 

Appendix  F  discusses  the  SMART  (  SMall  Acoustic  Ray  Trace,  Ref. 
22)  selected  as  the  eigenray  model  for  this  thesis*.  The  program  computes  the 
dependence  of  range  on  initial  angle  and  uses  this  relation  in  its  iterative  search 
for  the  eigenray. 

b.  Generation  of  Direct  Function  Table 

The  SMART  model  produces  the  sound  travel  times  for  an  eigenray 
set  pertaining  to  one  receiver  and  one  source.  The  generation  of  the  complete  direct 
function  table  requires  the  following  five  additional  steps. 

1.  Repeatedly  exercising  the  model  for  all  the  source  points  on  the  discrete  grid 
of  target  depth  and  range  (DR)  positions. 

2.  Selecting  the  shortest  travel  time  as  the  first  arrival  and  subtracting  it  from 
the  travel  time  of  the  other  paths.  The  resulting  differences  correspond  to  the 
time  lags  of  the  ACF. 

3.  Sorting  the  paths  in  ascending  travel  time  order.  This  corresponds  to  the 
estimator  logic  which  assumes  that  it  is  not  possible  to  directly  associate  the 
TDOAs  with  the  paths. 

4.  Applying  the  resolution  limitation  by  eliminating  delays  which  are  less  then 
the  resolvable  difference  away  from  their  preceding  delay. 

5.  Selecting  the  shortest  two  delays  and  storing  them  as  the  tabulated  direct 
function. 


The  last  step  is  in  need  of  further  elaboration.  The  data  structure  used 

for  the  purpose  of  storing  the  direct  function  table  is  a  3-D  array  Q  of  dimension 

2  x  Nd  x  Nr  indexed  by  £,  i,j.  The  first  array  index  selects  the  first  and  second 

*  The  proprietary  model  was  made  available  to  this  research  courtesy  of 
Mission  Sciences  Corp.,  Commack,  N.Y. 


delays  t\’,t2(t  =  1,2).  The  second  and  third  indices  t  and  j  select  a  point  on  the 
depth  and  range  grid  respectively.  The  variables  Nd  and  Nr  are  the  number  of 
depth  and  range  steps  in  the  grid.  Each  point  [t  j]  in  the  grid  thus  conceptually 
corresponds  to  a  four  element  vector  called  a  quad  ( q )  which  has  the  following  four 


values: 


Qij  ~  Dg  [i],  Rg  [i],  tl  [t,  j],  ^2  [*»  j] 


(3.22) 


where  Dg[i ]  is  the  depth  implied  by  the  depth  grid  index  t.  Rg\j]  is  the  range 
implied  by  the  range  grid  index  j.  t\,  t2  are  the  two  modified  TDOAs  given  by  tx  = 
T\—  To, and  <2  =  T2—T0,  and  stored  in  positions  Q[l,t,j]  and  Q[2,  i,j]  respectively, 
where  Tq,T\,  T2  are  the  ordered  travel  times  of  the  first  three  resolvable  paths  such 


that  T0  <  Ti  <  T2.  Namely 


t2  >  t j. 


(3.23) 


The  vectors  Dg  and  Rg  are  the  depth  and  range  grid  scale  vectors 


Dg  =  Dgo  +  AD  ■  (0, 1,2, . . .  Nd  -  1) 


Rg  =  Rgo  +  AR-  (0,1,2,... Nr  -  1 


(3.24a) 

(3.246) 


with  Dgo ,  Rgo  the  initial  grid  points  at  shallow  depth  and  short  range,  and  AD,  AR 
axe  the  grid  step  sizes. 

c.  Interpolation 

Linear  interpolation  of  the  TDOA  values  for  points  of  the  DR  grid  is 
used  to  provide  a  continues  TDOA  function  of  the  DR  position.  The  values  of  the 
*1(^2)  depth  and  range  of  the  three  quads 
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as  shown  graphically  in  Fig.  3.9. 
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Fig.  3.9.  Linear  interpolation. 


The  ordering  of  the  tabulated  discrete  direct  function  according  to 
depth  and  range  enables  finding  the  triangle  of  closest  quads  by  rounding  the 
input  D,nRtn  to  the  closest  DR  grid  values. 


d.  Implementation 


SMART  is  a  proprietary  product  which  was  obtained  for  this  research 
in  a  form  executable  on  an  IBM  personal  computer.  The  travel  time  data  generated 
by  this  program  was  transferred  to  a  large  mainframe  computer  (IBM).  Here,  the 
rest  of  the  processing  in  steps  2  through  5  and  the  target  tracking  simulations 
were  performed.  The  mainframe  was  selected  primarily  because  of  its  good  APL 
interpreter  and  the  excellent  engineering  graphics  support  provided  by  GRAFSTAT 
(Ref.  27]. 

3.  TDOA  of  IH  Surfaces 

Two  examples  of  the  resulting  interpolated  discrete  direct  function  surfaces 
are  presented  here.  The  form  is  the  same  as  that  used  for  the  homogeneous  cases 
presented  in  Section  B.  The  cases  are  coded  using  the  casename  code  described  in 
Appendix  E. 

An  inhomogeneous  finite  resolution  nonas9ociable  case  (A2251)  is  pre¬ 
sented  first.  The  TDOA  surfaces  of  SVP  with  a  positive  surface  gradient  of  0.05 
sec-1  and  very  high  delay  resolution  (0.05  ms)  are  shown  in  Fig.  3.10.  Note  that 
the  DR  grid  only  covers  the  first  100  m  of  the  500  m  water  column.  The  contour 
plot  of  the  1 1  surface  for  the  homogeneous  case  (case  S2251)  was  shown  in  Fig.  3.7. 

While  the  basic  trends  of  the  1 1  tj  surfaces  are  similar  for  homogeneous  and 
IH  cases,  the  specific  TDOA  values  differ.  This  can  be  seen  most  distinctively  on  the 
contour  plot  for  1 1  which  here  happens  to  be  the  delay  between  the  surface  bounce 
and  the  direct  path,  and  is  in  general  longer  than  the  t\  produced  by  straight  line 
propagation.  The  difference  results  from  the  fact  that  the  average  speed  of  sound 
is  smaller  along  the  IH  surface  bounce  path  than  along  the  homogeneous  surface 
bounce.  The  bottom  bounce  lag  #2  is  less  affected  since  the  most  of  its  path  is  in 
the  region  below  the  observer  where  the  speed  of  sound  is  constant  and  identical 


in  both  cases.  The  TDOA  surfaces  for  the  same  medium  conditions,  but  with  a 
finite  delay  resolution  of  0.5  ms  (Case  2251),  are  shown  in  Fig.  3.11.  The  inability 
to  resolve  the  surface  bounce  is  clearly  seen  along  the  first  depth  grid  line  of  both 
TDOA  plots. 

4.  2-D  Function  Inversion 
a.  Overview 

The  inversion  algorithm  is  an  on-line  transformation  of  a  2-D  TDOA 
input  measurement  vector  into  a  2-D  depth  and  range  position  vector  DR.  The 
time  delays  need  not,  and  in  most  cases  do  not ,  coincide  with  points  on  the  direct 
function  generating  grid. 

The  task  can  be  stated  graphically  when  the  2-D  transformations  are 
viewed  as  mappings  between  the  position  and  delay  planes.  The  direct  function 
maps  constant  depth  and  range  lines  (DR  grid)  into  curved  contours  on  the  tit2 
plane;  the  inverse  function  maps  a  txt2  grid  into  curved  contours  on  the  depth 
range  plane.  The  inversion  task  is  thus  defined  as  finding  the  intersection  of  the 
inverse  mapping  tx  t2  contours  and  reading  its  DR  coordinates. 

Typical  inverse  mappings  of  constant  <i<2  on  the  DR  plane  were  shown 
in  Fig.  3.12  graphically  demonstrating  the  inversion  process.  Two  contours  corre¬ 
sponding  to  t\  value  of  7.5  ms  and  <2  values  of  30  and  50  ms  are  shown  on  Fig. 
3.12.  One  intersection  of  the  7.5  and  30  ms  contours  is  shown  (a).  Four  intersection 
points  of  the  contours  corresponding  to  t\t2  values  of  7.5  ms  and  50  ms  jure  shown 
(b  through  e).  The  multiple  intersections  indicate  the  inherent  ambiguity  of  the 
realistic  case  which  was  discussed  in  Sec  B. 

Typical  direct  mappings  of  constant  depth  and  range  (DR  grid)  into 
the  <1*2  plane  axe  shown  in  Fig.  3.13a.  These  contours  are  computed  using 
ideal  straight  line  propagation.  From  Eq.  (2.106)  it  is  clear  that  the  depth  is 
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Fig.  3.11.  Inhomogeneous  and  unresolvable  TDOA  (C2251). 
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dependent  on  the  ratio  T\  /  r2  while  the  range  is  dependent  on  the  weighted  sum 
(Dw  —  D0)t\  4-  D0t2 .  This  dependence  explains  the  straight  line  nature  of  the 
inverse  direct  mapping.  The  dots  along  the  contours  of  Fig.  3.13  represent  quads. 
The  depth  and  range  are  defined  by  the  contours  and  the  TDOAs  are  the  Cartesian 
coordinates  of  the  dots. 

A  typical  direct  mapping  of  quads  for  am  IH  realistic  case  is  shown  in 
Fig.  3.13b.  Note  that  the  order  of  the  points  on  the  plane  is  lost  and  that  the 
points  are  restricted  to  the  half  plane  defined  by  the  <2  >  h-  This  results  from  the 
definition  of  the  modified  TDOA  (Recall  Section  B). 

A  cornerstone  of  the  inversion  is  the  fact  that  a  quad  qX}  is  in  itself 
invertible  in  the  sense  that  it  provides  the  transformation 


as  well  as  the  inverse  relation 


(3.31) 


(3.32) 


The  quads  can  therefore  be  viewed  in  reverse,  namely  as  defining  depth  and  range 
values  over  the  t\t2  delay  plane.  Linear  interpolation  can  be  called  upon  again 
to  compute  depth  and  range  values  at  f  j  t2  points  which  do  not  coincide  with  any 
quad. 

Although  the  inversion  procedure  seems  theoretically  quite  simple,  its 
practical  implementation  is  complicated  by  the  following  fact:  the  quads  in  the  Q 
array  are  not  ordered  on  the  tit2  plane  as  they  are  on  the  DR  plane  grid.  This 
is  especially  true  for  the  realistic  IH  case  as  shown  in  Fig.  3.13b.  The  search  for 
the  three  interpolation  quads  is  therefore  much  more  elaborate  than  the  simple 
rounding  procedure  used  for  the  direct  interpolation.  The  retrieval  of  the  quads 
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here  involves  a  search  which  is  drastically  effected  by  the  ambiguity  of  the  surfaces. 
The  search  is  described  in  the  next  section. 

Once  an  appropriate  quad  triangle  is  located,  depth  and  range  are 
interpolated.  If  multiple  solutions  are  found  one  has  to  be  chosen  and  reported 
as  an  output.  The  complete  procedure  includes  the  search,  the  interpolation,  and 
the  ambiguity  resolution.  The  inversion  is  defined  here  as  the  inversion  procedure 
(which  is  performed  by  the  prefilter).  The  inversion  procedure  is  described  in 
Section  C.4. 

b.  TDOA  Search 

Every  three  direct  function  quads,  each  relating  values  of  depth  D  (or 
range  R )  to  delays  t\t2  can  specify  a  plane  in  the  Dt\t2  (or  Rt\t2)  space  which 
can  be  used  for  the  inverse  interpolation.  In  order  to  improve  the  interpolation 
accuracy,  only  quads  forming  small  triangles  which  surround  the  input  t\t2  point 
are  considered.  The  total  number  of  combinations  of  three  quads  out  of  the  total 
(ND  x  NR)  quads  is  (jVD*  iVR)!  _  This  makes  exhaustive  search  for  the  smallest 
surrounding  triangle  an  impossibility. 

Since  different  triangles  produce  different  interpolation  results,  succes¬ 
sive  prefilter  iterations  for  the  same  tit2  point,  may  produce  different  depth  and 
range  values  if  the  selected  triangles  are  not  the  same.  This  is  an  undesirable  situ¬ 
ation.  For  example,  a  static  target  could  produce  changing  position  measurements, 
causing  the  target  tracker  to  respond  as  if  the  target  were  maneuvering. 

The  impact  of  a  nonunique  triangle  selection  is  demonstrated  in  Fig. 
3.14.  Part  c  and  e  of  the  figure  show  the  same  four  quads  (91929394)  grouped  in 
different  triangles  together  with  the  resulting  interpolated  range.  The  quads  are 
placed  on  a  cartesian  grid  for  clarity  purposes  only  and  the  values  marked  near  the 
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Fig.  3.14.  Non-unique  triangle  selection. 

quads  (8  and  12)  are  the  range  values  of  the  quads.  Three  interpolation  triangles 
can  be  selected  for  range  interpolation  for  a  target  moving  from  point  e  to  point  c. 
These  triangles  are  the  set  (case  A) 

Tril  :  (?3,  92,94) 


Tn2:  (<?i ,  <73 ,  <7i ) 


or  the  single  triangle  (case  B) 


TriZ  (91,73,94) 


Selecting  Tril  for  target  position  e  to  d  and  tn 2  for  position  d  to  c  (case  A)  will 
provide  a  continuous  range  interpolation.  Selecting  tn 3  for  the  complete  trajectory 
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a  to  c  will  also  yield  continuous  interpolation  (a  constant).  Selecting  Tn'3  for  posi¬ 
tions  e  to  d  followed  by  tri2  for  d  to  c  will  result  in  the  discontinuous  interpolated 
range  shown  in  d.  This  discontinuity  in  range  is  troublesome  for  the  maneuvering 
target  tracker. 

The  uniqueness  problem  has  two  aspects.  One  part  of  it  results  from 
the  unordered  nature  of  the  points  on  the  tit?  plane.  This  part  can  be  overcome 
by  constructing  unique  triangles  either  for  the  whole  grid  at  one  time  or  locally, 
when  and  where  needed  using  some  algorithm.  An  example  of  this  problem  is  the 
fact  that  every  triangle  that  surrounds  the  input  point  can  itself  be  surrounded  by 
another  yet  larger  triangle.  A  simple  rule  of  choosing  the  smallest  possible  triangles 
(for  example  by  area)  could  solve  this  problem. 

The  other  aspect  of  the  problem,  is  the  inherent  ambiguity  resulting 
from  the  non-monotonicity  of  the  direct  function.  This  ambiguity,  discussed  earlier, 
is  demonstrated  here  again  with  an  emphasis  on  its  relation  to  the  triangle  selection. 
Six  quads  marked  q\  through  q6  are  shown  mapped  on  the  t\t2  plane  in  Fig.  3.15a 
with  the  corresponding  depth  values  as  the  vertical  axis.  The  folding  surfaces 
formed  by  the  quads  helps  one  understand  the  name  Manifold  by  which  they  are 
celled.  As  expected,  the  non-monotonicity  of  the  direct  function  results  in  two 
possible  interpolation  triangles  frtl:  ( 91,92,93 )  and  tri2 :  (94,95,96)  for  an  input 
point-p.  The  problem  arises  because  both  triangles  are  “valid.”  Any  attempt  to 
resolve  the  dilemma  by  constructing  a  “unique”  triangular  grid  like  the  sample  in 
Fig.  3.15b  may  create  wrong  interpolating  triangles  by  mixing  delay  values  from 
two  separate  areas  of  the  non-monotonic  function  (e,g.  92,94,93). 

The  two  triangles  represent  two  separate  areas  on  the  DR  plane  for 
which  the  particular  MP  structure  generates  similar  delay  values  t\t2.  This  is 
caused  by  the  following  combined  effects:  the  IH  propagation;  the  limited  resolution 
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receiver;  and  primarily  the  lack  of  path  delay  association  and  the  corresponding 
assumed  instrumentation  logic  (described  in  Section  B).  The  triangles  selection 
mechanism  that  is  used,  has  to  retain  the  option  of  selecting  both  triangles  tril 
and  tri2  since  they  axe  both  valid  for  interpolation. 

Attention  is  directed  at  this  point  to  the  fact  that  while  tril  and 
tri2  are  partially  overlapping  in  the  TDOA  plane,  they  are  disjoint  in  the  DR 
plane.  This  difference  is  significant  and  is  used  later  to  assist  in  the  selection  of 
the  triangles. 

An  algorithmic  solution  to  the  general  problem  of  uniquely  selecting 
a  triangular  grid  over  a  set  of  unordered  points  in  a  plane  does  exist.  Surprisingly 
enough,  it  was  developed  only  recently  (1975)  by  Shamos  [Ref.  28].  The  solution  is 
elaborate  and  involves  a  geometric  construction  but  it  does  guarantee  unique  trian¬ 
gles.  However  it  does  not  solve  the  above  mentioned  problems  of  eliminating  valid 
and  creating  nonvalid  triangles.  The  Shamos  solution  is  therefore  not  applicable 
for  the  problem  of  non-monotonic  function  inversion  and  a  different  new  approach 
had  to  be  devised  here.  This  approach  is  described  next, 
c.  DR  Search 

The  difficulties  encountered  in  selecting  triangles  resulted  from  the 
unordered  scattering  of  the  quads  in  the  TDOA  plane  (Fig.  3.13b).  On  the  DR 
plane  the  quads  are  ordered  on  the  Cartesian  grid.  Constructing  a  unique  disjoint 
set  of  triangles  which  exhaustively  covers  the  depth  range  plane  is  trivial  and  can 
be  done  by  splitting  each  rectangle  of  the  DR  grid  into  two  triangles  using  parallel 
diagonal  lines  as  shown  on  Fig.  3.16.  The  set  of  triangles  thus  formed  is  defined 
as  the  DR  triangles  space.  A  key  idea  in  the  search  is  to  consider  only  triangles 
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Fig.  3.16.  Unique  valid  DR  triangle  set. 

in  the  delay  plane  that  arise  from  the  mapping  of  triangles  from  the  DR  triangles 
space.  This  delay  triangles  space  has  the  following  distinct  advantages: 

1.  The  set  is  unique. 

2.  The  set  is  completely  valid.  That  is,  the  triangles  Eire  locally  the  smallest  and 
will  not  be  constructed  by  erroneously  grouping  quads  from  disjoint  areas  of 
the  direct  function. 

3.  The  set  is  sufficient.  That  is,  the  option  of  choosing  all  valid  TDOA  triangles 
is  maintained. 

4.  The  search  space  size  is  significantly  reduced.  The  number  of  possible  triangles 

is  reduced  from  to  2  ■  Nd  •  Nr.  This  is  a  reduction  from  order  0(n!) 

to  order  O(n). 

5.  Construction  of  the  set  is  trivial  and  does  not  require  computation. 

6.  The  set  facilitates  an  efficient  structured  search.  The  search  for  the  surround¬ 
ing  TDOA  triangle  can  utilize  the  structure  of  the  underlying  DR  grid  in  an 
effective  binary  type  search. 


Viewed  in  the  context  of  the  graphical  statement  of  the  inversion  prob¬ 
lem  (Fig.  3.12)  the  search  now  is  for  the  DR  triangle  that  surrounds  the  intersection 
point  of  the  <!,£2  contours  of  the  inverse  function  on  the  DR  plane.  The  existence 
of  a  solution  is  conditioned  by  the  monotonicity  of  the  direct  function  over  the 
region  of  the  triangle  which  dictates  that  the  DR  grid  be  made  smaller  than  the 
minimum  monotonic  region. 

Since  the  DR  grid  of  quads  is  Cartesian,  searching  for  a  minimum  rect¬ 
angle  that  contains  the  desired  triangle  is  more  convenient.  After  the  surrounding 
rectangle  is  found  the  specific  triangle(s)  can  then  be  located  by  examining  the  two 
triangles  that  make  up  the  rectangle, 
d.  Localized  Search 

Finding  the  particular  triangle  in  the  total  set  of  2  -Nd-  Nr  grid  points 
that  surround  the  £j£2  input  measurement  point  requires  further  searching.  The 
test  to  see  if  the  point  is  within  a  triangle,*  involves  some  computation  and  the 
search  space,  though  reduced  can  be  still  very  large.  Hence  the  possibility  of  ex¬ 
haustively  considering  all  the  triangles  is  not  desirable.  A  binary-type  search  made 
possible  by  the  DR  ordered  grid  structure  can  reduce  the  average  computation 
load.  In  addition  it  can  focus  the  attention  on  the  area  where  the  solution  is  most 
likely  to  be  found. 

A  place  at  which  to  start  the  search  (a  root)  is  provided  by  the  target 
tracker  predicted  position.  If  a  target  state  vector  has  already  been  established, 
then  the  next  target  measurement  is  most  likely  to  be  found  in  the  vicinity  of  the 
predicted  DR  position.  Starting  the  search  with  the  rectangle  that  surrounds  the 
predicted  position  is  only  the  natural  choice. 


The  Q  array  indexes  required  to  retrieve  the  quads  associated  with 
this  rectangle  are  produced,  as  they  were  for  the  direct  function  interpolation, 
by  rounding  to  the  closest  DR  grid  values.  If  the  rectangle  that  surrounds  the 
predicted  point  also  surrounds  the  tx  t2  point  the  search  has  produced  the  required 
triangle  and  it  terminates.  Upon  failure,  however,  the  search  has  to  continue  over 
a  growing  and  more  distant  set  of  DR  rectangles  until  a  surrounding  one  is  located. 
This  local  search  ensures  that  the  solution  closest  to  the  target  will  be  found  first. 
This  significantly  reduces  the  ambiguity  inherent  in  the  MP  inversion. 

Two  different  search  strategies  were  designed  and  implemented  for  this 
purpose.  The  first,  titled  ME  A  (means  ends  analysis)  is  a  directed  search  strategy 
common  in  artificial  intelligence  applications.  The  second  method  referred  to  as 
local  search  algorithm  (LSA)  is  a  breadth  first  type  search  performed  in  tiers  of 
growing  distance  from  the  predicted  point. 

The  MEA  method  performed  very  effectively  for  monotonic  direct 
functions  but  failed  in  non-monotonic  cases.  The  LSA  algorithm  performed  very 
well  for  ail  cases  and  the  computational  load,  usually  associated  with  breadth  first 
search,  was  reduced  by  an  innovative  dual  phase  evaluation  algorithm.  A  brief 
description  of  the  local  search  algorithm  is  presented  next  and  a  more  detailed  ex¬ 
planation  of  its  components  is  presented  in  Appendix  G.  The  system  block  diagram 
is  shown  in  Fig.  3.17.  Note  the  feedback  of  the  predicted  position  as  a  root  for  the 
prefilter’s  local  search. 

e.  The  Local  Search  Algorithm  (LSA) 

Given  the  TDOA  measurement  point  the  goad  of  the  LSA  is  to  produce 
an  interpolation  triangle  that  surrounds  the  intersection  of  the  inverse  function  fif2 
contour  lines.  The  search  commences  in  a  DR  rectangle  where  the  solution  is  most 
likely  to  be  found,  namely  at  the  predicted  target  position.  The  search  progresses 


Fig.  3.17.  System  block  diagram. 

by  increasing  the  rectangular  search  area  until  a  solution  is  found  or  until  the 
dimensions  of  the  area  have  reached  a  specified  maximum  limit.  Upon  termination 
(success  or  failure)  the  algorithm  reports  the  number  of  solutions  found,  the  size 
of  the  area  actually  searched  and  the  triangles  found  (if  any). 

The  localized  search  is  a  part  of  the  overall  inversion  algorithm.  The 
inversion  initiates  and  uses  the  outputs  of  the  search  to  produce  the  desired  mea¬ 
surement  depth  and  range.  The  interface  between  the  inversion  and  the  search  is 
given  in  Table  3.1. 

A  numerical  analysis  procedure,  originally  developed  for  finding  zeros 
of  complex  functions  (Hamming  [29])  was  modified  to  apply  to  the  search.  One  test 
and  two  recursive  procedures  are  used  in  the  LSA,  they  are  the  inversion  test  the 
outward  and  the  inward  searches.  The  inversion  test  examines  a  given  rectangular 
on  the  DR  grid  to  find  out  if  a  solution  may  exists  in  it.  The  outward  search 
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TABLE  3.1 


LSA  INTERFACE  VARIABLES 

LSA  input: 

<1  <2  Delay  measurements  (input  point). 

DR  prd.  Target  predicted  position,  depth  range. 
Max  size  Maximum  search  area  size  , depth  range. 

LSA  output: 

Num  The  number  of  triangles  found. 

Area  size  The  size  of  area  searched  before 
the  solution  was  found. 

Solution  The  quad  triangles. 


examines  the  current  (input)  rectangle  using  the  inversion  test.  As  long  as  the  test 
fails  (no  solution  in  the  current  area)  and  as  long  as  the  maximum  search  area  has 
not  been  reached,  the  search  area  grows  outward.  If  the  maximum  area  is  reached 
and  a  solution  was  not  found  a  failure  is  reported. 

The  inward  search  attempts  to  locate  a  solution  inside  a  given  rect¬ 
angular  area.  It  starts  by  performing  the  conditional  inversion  test  on  the  entire 
area.  If  the  test  succeeds  the  area  is  divided  to  four  quadrants  and  a  solution  is 
searched  for  in  each  quadrant  by  means  of  a  recursive  call  to  the  inward  search. 
The  recursion  terminates  either  when  a  1  x  1  rectangle  is  reached,  or  when  the 
test  indicates  that  a  quadrant  does  not  contain  a  solution.  If  a  1  x  1  rectangle  is 
formed,  an  additional  surrounding  test  is  performed  on  the  two  DR  grid  triangles 
that  comprise  that  rectangle  to  determine  which  if  any  of  them  contains  a  solution. 
All  the  triangles  thus  found  are  reported  as  outputs  to  the  inversion  procedure  for 
further  processing. 

The  search  as  a  whole  is  thus  completed.  The  surrounding  triangles 
closest  to  the  predicted  target  position  are  reported  to  the  inversion  procedure 
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for  the  further  processing  which  is  discussed  next.  The  local  search  algorithm 
is  demonstrated  in  Fig.  3.18.  Fig.  3.18(a)  shows  two  iterations  of  the  outward 
search  initiated  at  point  p.  Fig.  3.18(b)  shows  three  iterations  of  the  inward  which 
produce  two  solutions  corresponding  to  the  two  intersections  of  the  t\tz  contour 
lines.  A  more  detailed  description  of  the  LSA  is  presented  in  Appendix  G  including 
a  detailed  pseudo  code  description, 
f.  The  Inversion 

The  inversion  procedure  performs  the  prefilter  operation  of  transform¬ 
ing  the  input  measurements  TDOAs  to  depth  and  range  measurements.  It  uses  the 
direct  function  table  and  the  local  search  algorithm  described  earlier.  The  inversion 
performs  the  following  five  tasks. 

•  Initialization  of  the  search 

•  Interpolation 

•  Ambiguity  resolution 

•  Failure  handling 

•  Performance  monitoring. 

Each  of  these  is  described  in  more  detail  below. 

(1)  Initiation  of  the  Search  Process.  The  process  is  initiated  by 
setting  up  the  first  lxl  search  rectangle  and  starting  the  outward  search.  The 
initial  rectangle  is  set  by  quantizing  the  target  predicted  position  to  the  closest  DR 
grid  position. 

(2)  Interpolation.  Once  the  appropriate  three  quads  are  located, 
the  depth  and  range  need  to  be  interpolated.  Linear  and  bilinear  interpolations 
were  used. 
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Linear  interpolation  is  done  in  much  the  same  way  as  the  time 
delay  interpolation  is  done.  The  interpolation  equations  axe  similar  to  those  de¬ 
scribed  earlier  for  the  direct  function  interpolation.  (See  Section  C.2.c.) 

A  linear  plane  in  the  depth  (range)  f  i ,  t2  space  of  the  form 

D  =  Adii  4-  Bdi2  +  Cd  (3.33) 

is  defined  by  the  three  quads  9*1  ,ji;  9*2.72;  9*3,73  where  9,^  is  defined  in  Eq.  (3.22) 
and  the  subscripts  index  the  three  points  of  the  selected  triangle  on  the  original 
depth  range  grid.  The  equations  for  the  plane  coefficients  in  the  depth  t\t2  space 
are 

■D[*i  ji]  —  Ad  -t\  [iiii]  4-  Bd  •  t2[ii  ji]  4-  Cd 

D[i2j2]  =  Ad  ■  ti  [i2j2]  4-  Bd  ■  t2[i2j2\  4-  Cd  (3.34) 

-0[i3 J3]  =  Ad  ■  *i  [13.73]  4-  Bd  •  t2 [13,73]  4-  Cd 
which  takes  the  matrix  form 

Ad 1  r<i[*iii]  l]  \Dlnh] 

Bd  =  *i  [12.72]  ^2  [*2  J2]  1  •  D[i2j2]  (3.35) 

.  Cd  _  tl  [*3^3]  ^2  [*3  J3]  1.  .Diizjz]. 

from  which  the  parameter  vector  ( AdBdCd)T  can  be  solved.  The  depth  (range)  at 
the  specific  t\mt2m  measurement  point  can  be  interpolated  as 

D  =  {AdBdCd)-(ti  t2  1)T.  (3.36a) 

And  similarly  for  range 

R  =  (Ar  Br  Cr)-(t  lm  t2rn  1  )r  (3.366) 

with  ArBrCr  solved  by  an  equation  like  Eq.  (3.35)  with  i2[i,j]  substituted  for 


Bilinear  interpolation  was  tested  as  well.  The  bilinear  interpo¬ 
lating  surface  [Ref.  30]  applied  to  depth  interpolation  takes  the  form 

D  =  Ad  •  <1  -f-  Bd  ■  <2  +  Cd  •  <i  •  t2  +  Dd-  (3.37) 

The  coefficients  Ad,  Bd,Cd,  Dd  of  the  bilinear  interpolation  can  be  determined  by 
solving  the  matrix  equation 

■<i[*iJi]  M*iii]  *1  [*1  ill  -  <2[*1  Ji]  1]  \Adl  rD[*,;j]' 

<1  [*2 J2 ]  <2 [*2 J2]  •  <2 [*272]  1  ,  Bd  _  Dfaji]  (V 

*  1  [* 3 J3 ]  Mz3j3  j  [*3j3]  *  <2[*3j3]  1  Cd  D[izjj] 

-*l[*4j4]  <2  [*4.74]  *l[*4i4]  -<21*474]  1J  \-Dd\  \-D\idji\. 

While  linear  interpolation  uses  a  triangle  formed  by  three  quads,  the  bilinear  inter¬ 
polation  requires  a  square  of  four  quads.  The  closest  square,  however,  is  naturally 
produced  by  the  local  search  algorithm. 

(3)  Ambiguity  Resolution  of  Multiple  Solutions.  When  multi¬ 
ple  solutions  exist  in  the  search  rectangle  a  choice  is  made  among  them  by  interpo¬ 
lating  all  of  the  candidate  solutions  and  selecting  the  one  closest  to  the  predicted 
position  (point  S2  in  Fig.  3.18).  For  a  single  linear  measurement  and  a  normally 
distributed  measurement  noise  and  prediction  error,  this  solution  will  be  the  more 
likely  one.  Realistically  the  measurement  is  nonlinear,  and  the  two  measurements 
(depth  and  range)  do  not  have  the  same  variance.  An  alternative  more  rigorous 
choice  would  imply  extensive  computations  that  are  expected  to  produce  practically 
the  same  results.  A  potentially  more  rigorous  approach  to  the  overall  ambiguity 
problem  is  discussed  in  Chapter  Six. 

(4)  Failure  Handling.  A  search  can  terminate  without  a  solution 
for  one  of  the  following  two  reasons: 

1.  The  measurement  noise  resulted  in  placing  the  ti  t2  input  point  outside  of  the 
maximum  search  area  or  even  outside  the  complete  DR  grid. 
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2.  The  search  area  reached  the  maximum  size  without  covering  the  actual  mea¬ 
sured  source  position. 

When  an  inversion  triangle  is  not  found,  the  predicted  target 
position  is  returned  by  the  inversion  procedure  as  an  output.  This  can  be  thought 
of  as  a  lack  of  innovation  in  the  current  measurement,  and  effectively  reduces  the 
sampling  rate  of  the  system.  The  tracker  performance  degrades  as  a  result.  When 
the  search  failure  occurs  infrequently  (low  failure  rate)  only  the  response  to  fast 
target  maneuvers  is  hampered.  Frequent  failures  (high  failure  rate)  can  lower  the 
effective  sampling  rate  of  the  observer  below  the  Nyquist  rate  required  to  sample  the 
target’s  dynamics.  This  drastically  reduces  the  quality  of  the  tracking.  Frequent 
failures  are  unfortunately  fed  back  in  the  form  of  search  areas  centered  around  the 
wrong  target  position.  This  in  turn  increases  the  chances  of  failure  and  eventually 
leads  to  a  complete  loss  of  tracking.  A  large  maximum  search  area  limit  will  reduce 
this  effect,  but  requires  increased  on-line  computation.  The  maximum  search  area 
limit  is  thus  a  key  system  parameter  which  needs  to  be  optimized  for  the  prevailing 
noise  conditions  and  the  available  computational  resources. 

In  experimental  simulation  of  the  MP  tracking  system  the  overall 
effect  of  the  inversion  failure  was  surprisingly  gradual  and  meaningful  tracking  of 
maneuvering  targets  was  demonstrated  with  up  to  a  50%  failure  rate.  The  robust 
performance  can  probably  be  attributed  to  a  more-or-less  symmetric  delay  error 
distribution.  In  such  cases  the  elimination  of  the  extreme  measurement  errors 
at  either  ends  of  the  distribution  caused  by  the  failure  does  not  affect  the  mean 
depth  and  range  values.  Therefore,  as  long  as  the  tracking  errors  are  not  large  and 
the  search  area  is  centered  properly  around  the  predicted  position,  the  impact  of 
eliminating  extreme  delay  measurements  is  mild. 
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(5)  Prefilter  Performance  Monitoring.  The  inversion  proce¬ 
dures  also  monitors  the  performance  of  the  search  by  averaging  the  number  of 
solutions  and  the  search  area  size.  The  averaging  is  done  by  means  of  an  auto¬ 
regressive  lowpass  filter.  The  average  solution  count  should  ideally  be  one,  a  lower 
value  indicates  search  failures,  while  a  higher  count  indicates  ambiguity.  The  search 
area  average  helps  optimize  the  maximum  search  area  and  is  a  measure  of  the  load 
associated  with  the  search.  It  also  provides  a  good  indication  of  the  effects  that 
delay  measurement  noise  has  on  the  search.  The  solution  count  and  the  area  size 
proved  to  be  very  useful  in  the  overall  system  performance  evaluation, 
g.  2-D  Function  Inversion  Summary 

A  2~D  function  inversion  algorithm  was  devised,  capable  of  translating 
realistically  estimated  MP  TDOA  generated  by  an  IH  medium  to  source  depth 
and  range  position.  An  eigenray  acoustic  model  identifies  and  precomputes  the 
the  travel  time  of  rays  to  the  receiver  from  hypothetical  targets  placed  at  the 
coordinates  on  a  depth  range  grid.  The  travel  time  is  transformed  to  TDOA  using 
an  assumed,  ACF  based,  realistic  receiver  and  delay  estimation  instrumentation 
and  stored  in  a  tabulated  direct  function.  The  inversion  is  performed  by  inverse 
interpolation  over  three  or  four  relevant  points  in  the  grid  located  by  means  of  a 
special  local  search  algorithm  (LSA).  A  total  of  three  prefilters  were  programmed. 
They  were  numbered  with  prefilter  numbers  (PFN)  4  through  6  (see  Chapter  Two, 
Section  E).  Prefilters  PFN-4  used  iineax  interpolation  and  the  MEA  search;  PFN-5 
used  linear  interpolation  and  the  LSA.  PFN-6  used  bilinear  interpolation  and  the 
LSA. 

The  procedure  was  carefully  evaluated  and  the  detailed  evaluation 
scheme  and  results  are  presented  in  the  next  chapter.  In  general  the  performance 


was  good  indicating  that  the  procedure 
path  tracking  scheme. 


An  intensive  evaluation  of  the  prefilter  was  performed  and  the  results  axe 
presented  in  this  chapter.  Section  A  discusses  the  motivation  for  the  evaluation 
and  some  of  the  expected  results.  Section  B  presents  the  evaluation  scheme.  Section 
C  contains  a  detailed  presentation  and  analysis  of  the  results, 

A.  MOTIVATION 

The  prefilter  forms  a  central  part  of  the  overall  MP  target  tracking  system. 
Examining  its  performance  is  a  necessary  prerequisite  to  analysis  of  the  overall 
system.  The  simulation  has  to  provide  all  the  realistic  conditions  under  which  the 
prefilter  is  expected  to  operate  as  a  part  of  the  MP  tracker.  In  addition,  the  sim¬ 
ulation  conditions  have  to  be  separately  controlled  in  order  to  isolate  their  impact 
on  the  prefilter.  The  effects  of  the  following  issues  on  the  prefilter  performance  are 
investigated: 

•  Medium  effects. 

o  IH  propagation, 
o  Propagation  loss, 
o  Bottom  depth  errors. 

•  Delay  estimation. 

o  Delay  noise  bias, 
o  Noise  range  dependence. 

•  Prefilter  design 

o  Grid  size. 


o  Interpolation  errors. 


The  overall  performance  of  the  prefilter  is  expected  to  vary  over  the  depth  and 
range  grid.  This  is  due  to  the  range  dependence  of  the  noise  and  the  nonlinearity 
of  both  the  inversion  and  the  time  delay  estimation  algorithms.  To  separate  the 
effects  the  evaluation  is  performed  over  the  entire  2-D  grid  in  two  steps.  In  the  first 
step  a  range  independent  TDOA  noise  is  assumed;  in  the  second  step  the  noise  is 
range  dependent. 

The  above  nonlinearities  give  rise  to  two  range  bias  errors.  The  first  is  a  range 
underestimation  bias  which  results  from  the  overestimated  TDOAs  produced  by 
1  *  nonnegative  effect  (see  Appendix  D.)  The  other  is  caused  by  the  delay  to 
position  transformation,  which  is  inherently  nonlinear.  This  second  bias  tends  to 
overestimate  the  range,  and  is  discussed  next.  The  effect  of  both  biases  on  depth 
is  dependent  on  the  particular  MP  structure. 

An  example  of  the  effect  of  the  nonlinear  TDOA  to  range  inverse  transforma¬ 
tion,  is  shown  in  Fig.  4.1.  A  section  of  the  range  dependence  on  tx  for  a  constant 
value  of  <2  is  shown  in  Fig.  4.1a.  A  zero  mean  normal  input  noise  distribution 
(Fig.  4.1c)  is  transformed  to  an  asymmetric  and  non  zero  mean  range  distribution 
(Fig.  4.1b).  The  typically  negative  value  of  leads  to  an  overestimated  range. 
This  can  be  seen  in  Fig.  4.1b  where  the  actual  range  is  1855  m  and  the  mean  is 
2050  m. 

Circles  representing  contours  of  a  joint  distribution  of  txt2  measurements* 
were  transformed  to  the  DR  grid.  The  results,  shown  in  Fig.  4. Id  demonstrate  the 
asymmetry  of  the  DR  measurement  noise. 

The  inversion  bias  was  studied  by  Moose  [11]  for  the  idealized  homogeneous 
case.  The  combined  opposing  effect  of  the  inversion  and  the  nonnegative  delay  are 
*  <7 =  a tj\tx  and  ^  independent. 


152 


PR0BA9U1Y  OtN'tnr  RANCt  [M) 


examined  here  for  the  IH  case.  Means  to  compensate  for  these  biases  are  discussed 
as  well. 


B.  THE  PREFILTER  EVALUATION  SCHEME 

The  evaluation  scheme  shown  schematically  in  Fig.  4.2  involved  the  following 
steps. 


Fig.  4.2.  Inversion  performance  evaluation. 

1.  Generation  of  the  reference  and  evaluated  direct  function  tables  using  the 
SMART  program  and  the  instrumentation  model  (see  Chapter  Three.).  Two 
separate  sets  of  conditions  are  used:  the  “actual”  (the  reference  case  or  Re- 
fcase)  and  the  “assumed”  (the  evaluated  case  or  Evalcase).  The  conditions 


such  as  SVP  bottom  depth  grid  size  etc.  may  be  different  for  the  Refcase  and 
the  Evalcase.  It  is  the  effect  of  these  differences  that  is  being  evaluated. 

2.  Generation  of  the  delay  measurement  noise  for  every  point  in  the  grid  and 
adding  it  to  the  delays  of  the  (reference)  TDOAs. 

3.  Simulation  of  the  noisy  target  predicted  position  Dp,  Rp  by  adding  depth  and 
range  noise  to  every  point  on  the  grid  values  Zfyj ,  R^j] .  This  supports  investi¬ 
gation  of  the  search  algorithm  in  an  open  loop  mode. 

4.  Transformation  of  the  measured  ti,<2  into  depth  and  range  measurement 
Dm,Rm  via  the  inversion  procedure  (prefilters  5  or  6)  using  the  evaluated 
direct  function  (i.e.,  the  assumed  ocean)  and  the  noisy  reference  point. 

5.  Comparison  of  the  measured  Dm,Rm  to  the  reference  Di^,Rij]  values  and 
generation  of  the  errors 

D error  =  Dm  -  D[i } 

Rerror  =  Rm  ^l\j\ 

6.  Repeating  the  above  for  all  the  points  in  the  DR  grid.  Recording  for  each  point 
the  errors  along  with  the  search  performance  measures  (number  of  solutions 
and  search  area  size). 

7.  Averaging  the  output  over  small  local  2-D  windows  to  provide  an  estimate  of 
the  local  mean  and  variance  and  a  local  prefilter  performance  measure. 

Recall  that  a  failed  search  will  return  the  reference  point  as  the  inverted  out¬ 
put.  When  noise  is  not  added  to  the  reference  point  this  will  produce  a  zero  error 
which  is  not  a  precise  account  for  the  event.  This  will  become  evident  from  some 
of  the  results. 

C.  EVALUATION  RESULTS 

The  evaluation  was  performed  about  150  times  for  a  variety  of  conditions.  The 
results  of  eight  summary  cases  are  presented  next.  The  cases  investigate  the  effects 
of  the  following  issues: 


•  Depth  range  grid  size. 

•  Linear  versus  bilinear  interpolation. 

•  TDOA  bias  error. 

•  Acoustic  medium  inhomogenuity. 

•  Range  independent  TDOA  noise. 

•  Range  dependent  TDOA  noise. 

•  Range  dependent  prediction  point  reference. 

•  Bottom  depth  error. 

Each  of  the  next  eight  subsections  describes  one  evaluation.  Each  subsection 
describes  the  purpose  of  the  simulation  experiment  and  the  expected  results,  the 
actual  experimental  results  and  analysis  of  the  evaluation,  and  conclusions.  The 
output  plots  include: 

•  Local  average  of  depth  and  range  errors. 

•  Local  average  of  inversion  error  STD. 

•  Local  average  of  the  solution  count  (this  monitors  the  performance  of  the 
search  algorithm.) 

•  Average  error  as  function  of  range  (range  error  profile).  This  is  generated  by 
taking  error  averages  over  all  of  the  depth  points  for  a  given  range.  The  error 
profile  is  normalized  as  marked  on  the  plots. 

1.  Effects  of  Grid  Size 

The  purpose  of  this  Run  was  to  examine  the  effect  of  grid  size  on  interpo¬ 
lation  error  and  optimize  the  grid  for  the  intended  use.  The  key  parameters  used 
are  given  in  Table  4.1. 

The  depth  and  range  results  of  the  inversion  with  linear  interpolation  over 
the  medium  density  DR  Evalcase  grid  (Si 111),  are  compared  to  the  exact  depth 
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TABLE  4.1 


INTERPOLATION  EVALUATION  GRIDS 


Grid  density 

High 

Medium 

Low 

Casename 

S1155 

Sllll 

S1144 

Initial  depth  [mj 

160 

160 

160 

Depth 

Step  size  [m] 

10 

50 

100 

Number  of  steps 

11 

11 

6 

Initial  range  [m] 

500 

500 

1000 

Range 

step  size  [m] 

250 

500 

1000 

number  of  steps 

99 

25 

13 

and  range  values  of  the  high  density  DR  Refcase  (Si  155)  grid. 

The  results  are 

shown  in  Fig.  4.3.  The  maximum  errors  were  16  m  in  depth  and  100  m  in  range. 

An  interesting  point  was  revealed  by  these  plots.  The  interpolation  er¬ 
rors  depend  on  the  degree  of  curvature  in  the  direct  and  inverse  functions.  The 
curvature  is  very  small  at  long  ranges  making  the  interpolation  error  inversely 
proportional  to  the  range.  This  is  clearly  seen  in  Fig.  4.3.  The  significance  of 
the  interpolation  errors  is  greater  at  short  ranges  than  at  long  ranges  since  the 
other  inherent  MP  tracking  errors  axe  reduced  at  short  range.  The  grid  size  should 
therefore  be  optimized  for  the  short  range  performance. 

A  more  coarse  grid  (Evalcase  Si  144)  was  evaluated  using  the  medium 
density  grid  as  the  reference  (Refcase  Sill).  The  results  shown  in  Fig.  4.4  reach  a 
maximum  depth  error  of  20  m  and  a  maximum  range  error  of  200  m. 

Conclusion:  The  grid  size  depends  primarily  on  the  direct  function.  A 
more  curved  and  more  rapidly  changing  function  will  require  a  finer  grid  to  maintain 
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low  errors.  For  the  nominal  homogeneous  case,  in  water  depth  on  the  order  of  500 
m  and  ranges  up  to  about  20  km  a  direct  function  Q  grid  size  of  2  x  20  x  40  will 
be  required.  This  is  a  very  reasonable  size  which  is  easily  within  the  capability  of 
current  onboard  microcomputers. 

2.  Bilinear  Interpolation 

The  purpose  of  this  evaluation  was  to  compare  the  performance  of  linear 
and  bilinear  interpolation.  The  key  parameters  used  are  the  same  as  those  used  for 
the  medium  grid  evaluation  (See  Table  4.1)  with  S1155  as  the  Refcase  and  Sllll 
as  the  Evalcase. 

Error  surface  for  the  bilinear  interpolation  are  plotted  in  Fig.  4.5.  The 
difference  between  the  results  of  the  two  interpolation  methods  is  small  (compare 
to  Fig.  4.3).  The  smaller  depth  errors  in  the  linear  interpolation  result  from  the 
fact  that  the  triangular  area  of  interpolation  is  smaller  than  that  of  the  rectangular 
area  used  by  the  bilinear  method.  The  bilinear  interpolation  provides  a  smoother 
interpolation  but  has  a  larger  absolute  error. 

Conclusion:  Both  linear  and  bilinear  interpolation  produce  satisfactory 
results.  The  error  of  the  linear  method  is  slightly  smaller.  The  choice  of  which 
method  to  use  cm  be  based  on  implementation  considerations.  Linear  interpolation 
was  used  for  most  of  this  study. 

3.  The  Effects  of  Acoustic  Medium  Inhomogenuitv 

The  purpose  of  this  evaluation  run  was  to  examine  the  effects  of  the 
medium  inhomogenuity.  A  SVP  with  a  positive  surface  gradient  of  steepness  that 
is  characteristic  of  many  practical  problems  ( g  =  0.05  sec-1.)  was  used  as  the  ref¬ 
erence  case  (A2251)  and  the  TDOAs  were  converted  to  depth  and  range  assuming 
a  homogeneous  medium  (straight  line  propagation,  Evalcase  S2251). 


Depth  and  range  inversion  errors  axe  shown  in  Fig.  4.6  as  surfaces  and  as 
contour  plots.  Very  large  errors  are  developed  especially  at  long  ranges.  Observe 


that  at  a  range  of  12  km  the  range  error  is  4  km  and  the  depth  error  is  150  m. 

The  trend  of  the  error  is  expected  and  is  caused  primarily  by  the  change  in 
the  surface  bounce  path.  Fig.  4.7  shows  a  comparison  of  an  IH  and  a  straight  line 
propagation  where  and  t2  for  both  cases  are  the  same.  The  path  T\  is  longer 
(curved  instead  of  straight)  and  includes  portions  with  average  speed  of  sound 
lower  than  in  the  homogeneous  case.  The  resultant  delay  t\  given  by  T\  —  To,  is 
longer. 

When  interpreted  as  if  it  were  the  result  of  straight  line  propagation  it  gives 
rise  to  deeper  depth  and  shorter  range  estimates. 

Conclusion:  The  inhomogenuity  of  the  ocean  medium  has  a  gross  effect 
on  MP  propagation  and  position  measurement.  Compensating  for  this  effect  is 
mandatory  in  order  to  achieve  meaningful  depth  and  range  measurements. 

4.  Finite  Delay  Resolution 

This  case  examined  the  effect  of  finite  delay  resolution  combined  with  the 
effect  of  the  inhomogenuity.  The  cases  used  were  Ref  case  C4261  for  which  the 
assumed  delay  resolution  is  0.5  ms  and  Evalcase  S4261  for  which  the  assumed 
delay  resolution  is  perfect  and  the  assumed  SVP  is  homogeneous. 

The  errors  are  similar  in  nature  to  the  errors  produced  by  the  previous 
run.  The  main  difference  is  along  the  surface  as  -hown  on  Fig.  4.8.  This  difference 
results  from  the  fact  that  in  the  finite  resolution  case  the  surface  bounce  is  not 
resolvable  along  the  first  depth  line  of  the  DR  grid.  The  first  depth  grid  point  is 
located  at  a  depth  of  2m  from  the  surface*.  The  delay  between  the  direct  path  and 

*  The  observer  depth  is  162  m  and  the  grid  initial  depth  is  +160  m  relative  to 
the  observer. 
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the  first  surface  bounce  is  smaller  than  0.5  ms  and  therefore  is  not  resolvable.  As  a 
result  the  modified  TDOA  vector  t\ ,  is  comprised  of  ri ,  r3  here  (instead  of  T\ ,  r2 ). 
This  corresponds  to  a  physically  different  multipath  configuration*.  Similar  effects 
were  demonstrated  near  the  ocean  bottom  and  along  the  DR  loci  where  t\  =  t2. 

Note  also  that  the  magnitude  of  the  error  over  the  complete  DR  grid  is 
larger  than  in  Fig.  4.6.  This  results  from  the  steeper  surface  gradient  of  0.105 
[sec-1]  (SVP  4  in  Appendix-E)  used  here  compared  to  the  0.05  [sec-1]  gradient 
used  in  Fig.  4.6. 

Conclusion:  The  finite  delay  resolution  has  a  significant  effect  on  the 
TDOA  structure.  Accounting  for  the  specific  response  of  the  delay  estimator  to 
the  loss  of  resolution  is  mandatory,  especially  for  targets  near  the  ocean  boundaries. 


The  purpose  of  this  evaluation  is  to  examine  the  effect  of  a  constant  TDOA 
offset.  This  was  done  as  a  scaling  run  and  supported  the  independent  evaluation  of 
the  tracker.  It  allowed  nominal  transformation  of  the  predicted  delay  noise  to  depth 
and  range  measurement  errors.  Since  this  is  a  nominal  scaling,  the  homogeneous 
case  (S3252)  is  used  both  as  the  Refcase  and  the  Evalcase.  A  ^  offset  of  1  ms  is 
used  as  TDOA  bias. 

The  resulting  error  is  shown  on  Fig.  4.9  and  reaches  levels  of  up  to  4  km 
in  range  and  80  m  in  depth  at  range  of  25  km.  With  actual  range  in  the  vicinity 
of  12  km  the  range  error  is  1  km  and  the  depth  error  is  40  m. 

Conclusion:  If  the  tracker  operates  over  a  range  of  up  to  12  km  and 
TDOA  noise  of  the  order  of  1  ms  is  assumed,  the  range  error  STD  can  be  expected 
to  lie  between  a  few  hundred  meters  and  about  1  km.  As  was  shown  in  Chapter 
Two,  such  errors  cm  be  handled  by  the  MM  target  tracker. 

*  See  Chapter  Three,  Section  III.B.5. 


6.  Effects  of  Range  Independent  Noise 

The  effect  of  a  stationary  (range  independent)  random  noise  on  both  the 
variance  and  the  mean  of  the  DR  measurement  was  examined.  The  Refcase  and 
Evalcase  were  both  taken  as  the  homogeneous  S4261.  A  stationary,  range  indepen¬ 
dent  (( p  =  0),  see  Eq.  (3.3))  TDOA  noise  with  STD  of  1  ms  for  both  t\  and  £2  was 
used.  The  local  averaging  window  for  the  plots  was  3  x  7  in  grid  step  units. 

Surface  plots  of  locally  averaged  STD  of  the  depth  and  range  errors  are 
shown  on  Fig.  4.10a.  The  depth  and  range  errors  were  averaged  over  depth  to 
provide  the  range  error  profile  which  is  shown  on  Fig.  4.10.  The  profiles  are 
normalized  to  IOC  m  in  depth  and  1  km  in  range. 

The  maximum  error  is  consistent  with  the  earlier  scaling  results  performed 
using  the  constant  TDOA  offset  (Fig.  4.9).  Contours  of  locally  averaged  number 
of  solutions  (Fig.  4.10d)  indicate  frequent  failure  of  the  inversion  near  the  deep 
end  of  the  grid*.  This  is  caused  by  the  TDOA  values  produced  from  the  noisy 
measurements.  The  noise  drives  the  TDOA  measurements  out  of  the  interval  of 
TDOA  values  covered  by  the  precomputed  direct  function.  Increasing  the  area 
covered  by  the  DR  grid  will  solve  part  of  the  problem  and  reduce  the  frequency  of 
failures. 

The  results  of  a  similar  case  with  TDOA  error  STD  of  5  ms  axe  shown  in 
Fig.  4.11.  This  figure  shows  larger  errors  and  failure  of  the  algorithm  at  a  shorter 
range  (7  km  in  Fig.  4. lid).  Note  that  the  term  failure  indicates  only  that  the 
effective  sampling  rate  is  reduced  by  half  at  ranges  over  7  km  (recall  the  comment 
at  the  end  of  Section  B).  More  frequent  failure  is  the  reason  for  the  superficially 


*  The  frequency  of  failure  is  50  %  along  the  0  solution  contour  line  since  the 
line  separates  the  regions  between  0  and  1  average  solutions.  The  contour  line  does 
not  imply  that  there  axe  no  solutions  at  ranges  longer  than  those  indicated  by  the 
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lower  error  values  at  the  long  range.  The  effect  of  the  reduced  sampling  rate  on 
the  overall  system  is  investigated  in  Chapter  Five. 

Conclusions:  Errors  with  STD  on  the  order  of  30  m  in  depth  and  0.5  km 
in  range  are  expected  at  ranges  around  10  km  as  a  result  of  a  range  independent 
TDOA  noise  with  STD  of  1  ms.  For  TDOA  noise  with  STD  of  5  ms  similar  depth 
and  range  errors  are  expected  at  a  range  of  5  km.  As  shown  in  Chapter  Two  these 
errors  can  be  handled  by  the  target  tracker.  The  inversion  algorithm  does  fail  at 
long  ranges  if  the  noisy  delay  measurements  fall  outside  of  the  interval  covered  by 
the  precomputed  direct  function. 

7.  Effects  of  Ranee  Dependent  Noise 


The  response  of  the  inversion  of  TDOA  which  is  contaminated  by  range 
dependent  noise  was  examined  here.  Range  dependent  noise  was  added  to  the 
target  predicted  point  as  well.  Case  S4261  was  used  both  as  the  Refcase  and  the 
Evalcase.  Range  dependent  noise  (Eq.  (2.96)  NT  =  1,  P  =  1),  with  ti,t2  STD  of 
0.2  ms  at  initial  range  of  500m,  was  used.  The  TDOA  noise  STD  varied  linearly 
with  ranges  between  0.2  ms  to  4.8  ms  over  the  range  interval  of  0.5  to  12.5  km. 
This  noise  level  resembles  the  one  produced  by  the  more  elaborate  noise  model 
(NT=2)  for  a  SNR  of  60dB  and  p  =  2  which  is  a  fairly  typical  case.  Predicted 
position  DR  noise  was  varied  from  0.4  to  250  m  in  depth  and  from  10  to  6.25  m  in 
range  proportional  to  the  second  power  of  range  (pj  =  2). 

The  resulting  error  profile  is  shown  on  Fig.  4.12  along  with  the  contour 
line  for  the  number  of  solutions.  Note  the  sharp  range  dependence  of  the  inversion 
error  STD  on  range.  Also  note  that  for  the  practical  case  used  here  a  sharp  “knee” 
develops  around  12  km  (Fig.  4.12c).  Beyond  this  range  the  error  grows  very 
rapidly.  This  close  to  second  order  dependence  of  the  DR  measurement  noise  on 
range  is  what  motivated  setting  the  parameter  pi,  to  2.  The  maximum  area  size 
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Fig.  4.11.  Range  independent  TDOA  noise  (5  ms). 


was  set  to  20  x  25  grid  points.  This  explains  why  the  addition  of  the  DR  noise  did 
not  degrade  the  performance  significantly.  Performance  degradation  as  a  result  of 
smaller  maximum  area  size  is  demonstrated  in  Chapter  Five. 

Note  that  the  normalized  range  STD  in  Fig.  4.12b  at  a  range  of  9  km  is 
approximately  0.3.  This  implies  a  range  error  STD  of  3.6  km  (considering  the  50% 
failure  rate  and  the  6  km  normalizing  factor).  This  value  will  be  useful  in  Chapter 
Five  in  setting  the  maximum  area  size. 

A  second  interesting  observation  is  the  failure  along  the  first  range  grid 
line.  This  is  caused  by  the  noise  which  drives  the  “measured”  TDOA  outside  of 
the  interval  covered  by  the  precomputed  direct  function.  Recall  however  that  the  0 
contour  line  in  the  number  of  solutions  plot  indicates  only  that  the  average  failure 
rate  is  larger  than  50%.  Similar  phenomena  were  noticed  in  most  simulations  along 
the  edges  of  the  DR  grid. 

Conclusions:  The  inversion  varies  with  range  even  for  a  range  indepen¬ 
dent  TDOA  noise.  The  error  becomes  even  more  range  dependent  when  the  TDOA 
noise  itself  is  range  dependent.  Beyond  15  km  very  large  DR  measurements  develop 
even  for  moderate  initial  range  delay  estimation  noises. 

8.  Effect  of  Bottom  Depth  Error 

The  purpose  of  this  last  inversion  evaluation  was  to  quantitatively  examine 
the  effects  of  an  erroneous  bottom  depth  on  the  target  depth  and  range  errors. 
Refcase  S1256  for  which  *He  bottom  depth  is  513  m  was  compared  to  Evalcase 
S1556  for  which  the  bottom  depth  is  550  m.  A  homogeneous  ocean,  perfect  delay 
resolution  and  no  noise  were  used  here  in  order  to  isolate  the  effects  of  the  erroneous 
bottom  depth.  The  inversion  was  done,  however,  using  the  numerical  algorithm 
(prefilter- 5). 
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The  effect  of  this  relatively  large  bottom  depth  error  (37  m  which  is  7.2% 
of  the  total  ocean  depth)  is  shown  on  Fig.  4.13  for  short  and  long  ranges.  While 
the  resulting  range  error  reached  a  level  of  3  km  at  18km,  it  was  limited  to  1.5  km 
for  ranges  smaller  then  10  km. 

The  trend  of  the  error  surfaces  is  explained  as  follows.  The  actual  (Refcase) 
ocean  produces  a  bottom  bounce  TDOA  which  is  shorter  then  the  one  produced  for 
the  same  DR  point  by  the  assumed  (Evalcase)  ocean.  This  smaller  ti  is  interpreted 
by  the  inversion  (  which  is  done  using  the  Evalcase)  as  the  TDOA  of  a  deeper  target 
at  a  longer  range.  The  region  of  no  inversion  solution  at  ranges  in  excess  of  18  km 
resulted  from  the  Refcase  <2  TDOA  being  smaller  then  the  minimum  #2  in  the 
precomputed  Evalcase  Table. 

Conclusions:  Incorrect  bottom  depth  assumptions  can  give  rise  to  large 
errors  in  MP  position  measurements,  primarily  at  long  range.  The  trend  of  the 
error  is  consistent  with  other  MP  tracking  error  sources  in  that  it  increases  with 
increasing  range.  This  error  tends  to  limit  practical  tracking  in  shallow  water  to 
ranges  below  20km. 
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Fig.  4.13.  Effort  of  bottom  depth  error. 


D.  SUMMARY  OF  PREFILTER  EVALUATION 


The  results  of  these  experiments  can  be  summarized  as  follows. 


1.  The  effect  of  the  inhomogenuity  on  the  MP  position  measurement  is  severe 
and  was  clearly  demonstrated. 

2.  The  depth  and  range  measurement  errors  are  strongly  dependent  on  range 
(proportional  to  2nd  to  4tA  power  of  range  depending  on  SNR  and  propagation 
loss). 

3.  Errors  with  STD  of  up  to  40  m  in  depth  and  1  km  in  range  can  be  expected 
for  ranges  of  about  7  km  and  typical  targets  with  SNRo  =  60  dB  if  spherical 
spreading  is  assumed  ( p  =  2).  Above  12  km  the  errors  increase  sharply  with 
rang~  f'  -  the  typical  shallow  water  case  studied  here. 

4.  The  inversion  algorithm  fails  under  very  noisy  conditions  and  around  the  edges 
of  the  DR  grid.  The  failures  in  these  cases  are  legitimate  since  the  TDOA  input 
values  are  not  representative  of  points  in  the  DR  grid. 

5.  Errors  will  result  if  account  is  not  taken  for  the  loss  of  delay  resolution.  In 
particular  one  needs  to  model  the  response  that  is  expected  from  the  specific 
type  of  delay  estimator. 

6.  The  algorithm  is  sensitive  to  DR  grid  step  size  only  at  short  ranges.  In  general 
the  depth  and  range  step  sizes  should  be  tailored  to  the  fluctuation  of  the 
specific  direct  function  used  and  the  accuracy  desired  at  short  range.  Grids  of 
10  to  50  m  in  depth  and  250  to  500  m  in  range  were  used  in  most  of  our  work. 

7.  The  MP  tracking  is  sensitive  to  bottom  depth  errors  primarily  at  long  ranges. 
At  ranges  below  10  km  a  low  percentage  error  in  bottom  depth  would  produce 
about  twice  that  percentage  error  in  range. 


V.  OVERALL  PERFORMANCE  EVALUATION 


A.  INTRODUCTION 

In  this  chapter  the  overall  performance  of  the  MP  tracking  system  is  evalu¬ 
ated.  A  total  of  seven  simulation  runs  which  axe  divided  into  into  three  groups 
are  presented.  The  runs  are  numbered  6  through  12  following  the  sequence  of  the 
runs  begun  in  Chapter  Two.  Runs  6  and  7,  represented  in  section  B,  demonstrate 
the  functional  performance  and  limitations  of  the  integrated  prefilter  and  tracker. 
Runs  8  through  10,  described  in  Section  C,  were  designed  to  investigate  the  effects 
of  medium  inhomogenuity  and  the  ability  of  the  new  prefilter  to  compensate  for 
it.  Runs  11  and  12,  covered  in  Section  D,  examine  the  sensitivity  of  the  IH  com¬ 
pensation  to  the  accuracy  of  SVP  measurement.  Conclusions  are  summarized  in 
Section  E. 

The  results  axe  presented  in  a  format  similar  to  the  one  used  for  the  previous 
runs  and  case  evaluations  in  Chapters  Two  and  Four.  In  each  section  the  goal 
of  the  runs  and  the  key  parameters  used  are  described  first.  Then  the  analysis 
results  are  presented.  The  detailed  parameter  values  for  all  of  the  runs  is  given  in 
Appendix  E. 


B.  INTEGRATED  SYSTEM  FUNCTIONAL  PERFORMANCE 

The  overall  performance  of  the  prefilter  integrated  in  a  closed  loop  with  the 
target  tracker  was  demonstrated  using  a  maneuvering  target.  The  target  initially 
has  an  outgoing  speed  (range  rate)  of  6  m/sec  and  it  makes  a  sharp  maneuver  at  a 
range  of  9  km  turning  back  towards  the  receiver  at  the  same  speed.  Realistic  noise 
(noise  type  2,  SNR  =  50  dB  and  p  =  2)  is  used  to  investigate  the  overall 


capability  of  the  system  through  the  maneuver.  A  constant  speed  of  sound  is  used 
both  as  a  reference  and  evaluation  medium  (case  Si  1 1 1 ).  This  is  done  in  order 
to  concentrate  on  the  functional  performance  of  the  integrated  system  when  the 
measurements  are  noisy. 

The  key  issues  investigated  here  relate  to  the  inversion  method.  These  issues 
are  the  ability  of  the  local  search  algorithm  (LSA)  to  find  a  valid  surrounding 
triangle  and  the  effects  that  a  reduced  effective  sampling  rate  has  on  the  target 
tracker  (recall  the  discussion  of  inversion  failure  in  Chapter  Three). 

Results  of  Run  6  are  shown  in  Fig.  5.1.  The  maximum  search  area  size  was 
±  200  m  in  depth  and  ±3.75  km  in  range  around  the  predicted  target  position 
(a  rectangle  of  8  x  15  DR  grid  units  for  case  S 1 1 1 1 ).  Fig.  5.1  shows  that  the 
overall  tracking  performance  is  good.  The  range  tracking  error  reduces  to  around 
300  m  at  the  14th  minute  after  initialization  is  completed.  The  realistic  TDOA 
noise  simulated  here  increases  rapidly  with  range.  This  is  clearly  evident  in  Fig. 
5.1  (note  the  two  resulting  range  measurements  errors  around  the  13tA  minute). 

The  system  tracks  well  through  the  maneuver  with  an  effective  sample  rate  as 
low  as  40%.  This  is  indicated  by  the  average  number  of  solutions  which  is  close 
to  0.4  in  Fig.  5. Id.  Note  that  the  prefilter  uses  the  maximum  search  area  size 
during  the  17th  minute  as  indicated  by  the  normalized  size  in  Fig.  5. Id  (actual 
searched  area  /  maximum  area  size).  Some  ambiguity  is  indicated  at  shorter  ranges 
where  the  soh’tion  count  is  larger  than  1.  This  ambiguity  affects  mostly  the  depth 
tracking  as  can  be  seen  in  Fig.  5.1b.  This  effect  is  expected  when  one  reviews 
the  TDOA  contours  of  case  Sllll  which  is  used  here  (see  Fig.  3.6).  The  bearing 
channel  tracking  and  the  resultant  XY  plot  of  Run  6  are  shown  for  completeness 
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Fig.  5.2.  Run-6:  Bearing  and  XY  plot. 

The  same  scenario  was  used  in  Run  7  with  a  smaller  maximum  search  area  size 
of  ±125  m  in  depth  and  ±2.25  km  in  range  (a  rectangle  of  5X9  in  the  Si  111  grid). 
The  results,  shown  in  Fig.  5.3,  clearly  indicate  the  failure  of  the  system  to  track  the 
target  through  the  maneuver.  The  system  loses  track  when  the  range  tracking  error 
reaches  2  km  (in  the  20* ^  minute).  At  this  time  and  range  most  of  the  measurements 
are  outside  of  the  maximum  search  area  around  the  predicted  position.  This  leads 
to  a  constant  use  of  the  maximum  search  area  size  and  eventually  to  a  complete 
failure  (over  80%  of  the  iterations). 

The  range  measurement  error  STD  for  Runs  6  and  7,  was  estimated  using  the 
inversion  performance  evaluation  scheme  described  in  Chapter  Three.  The  typical 
STD  for  similar  TDOA  noise  at  around  9  km  was  3.6  km*.  The  conclusion  is 
*  Recall  the  discussion  in  Chapter  Four,  Section  C.7 


therefore  that  the  search  area  size  should  be  set  to  twice  the  standard  deviation  of 
the  measurement  noise,  ((2crj)  x  (2ay)). 

C.  TRACKING  IN  AN  INHOMOGENEOUS  MEDIUM 

The  next  three  Runs  demonstrate  the  ability  of  the  new  MP  tracking  scheme 
to  compensate  for  the  effects  of  ocean  inhomogenuity.  A  reference  case  (Refcase) 
was  used  to  generate  the  TDOAs  (the  direct  function  table)  used  for  simulation  of 
the  “actual”  medium.  The  precomputed  direct  function  table  generated  using  the 
evaluated  case  (Evalcase)  was  used  in  the  prefilter.  The  Evalcase  represented  the 
-ssumed  ocean  conditions. 

Rim  8  demonstrates  the  effects  of  an  IH  medium  (Refcase  C2251)  when  straight 
line  propagation  is  assumed  and  used  in  the  prefilter  (Evalcase  Sllll)  .  The  results 
are  shown  in  Fig.  5.4.  The  large  range  tracking  error  of  3.3  km  at  a  range  of  10  km 
clearly  indicates  the  need  for  compensation  for  the  medium  inhomogenuity.  A 
relatively  strong  target  signal  (SNRo  70  dB)  was  used  in  order  to  isolate  the  IH 
medium  effects  from  the  effects  of  the  noise. 

Run  9  examines  the  same  scenario  with  a  weaker  target  signal  (SNRo  =  50 
dB)  and  its  results  are  shown  in  Fig.  5.5.  The  TDOA  measurements  are  noisy  now 
and  give  rise  to  a  negative  range  bias  which  adds  to  the  effect  of  the  IH  ocean.  The 
total  resulting  estimation  error  for  this  case  is  4  km  at  the  10  km  range. 

The  use  of  the  new  prefilter  with  the  correctly  assumed  IH  medium  (Evalcase 
C2251)  is  demonstrated  in  Run  10  and  the  results  are  shown  in  Fig  5.6.  The 
removal  of  the  error  associated  with  the  IH  medium  is  clearly  evident.  The  error 
in  the  depth  channel  after  the  depth  maneuver  at  the  8th  minute  is  caused  by 
the  fact  that  the  C2251  grid  was  limited  to  a  depth  of  -f  60m  (due  to  the  current 
limitation  of  the  SMART  eigenray  program  discussed  in  Appendix  F).  This  error 
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Fig.  5.4.  Run-8:  Tracking  in  an  IH  medium. 

is  not  present  at  shallower  depths  away  from  the  grid  edge  or  if  the  noise  is  smaller 
as  demonstrated  by  the  near  perfect  tracking  in  Fig.  5.6  where  SNR0  of  70  dB 
was  used.  This  error  can  be  removed  completely  if  a  full  grid  is  employed.  The 
average  number  of  solutions  reached  1.8  in  some  cases  since  the  corresponding 
TDOA  contour  lines  intersect  in  four  closely  located  DR  positions  (see  Fig.  3.17 
which  is  taken  from  typical  contours  of  the  same  case  C2251).  The  effects  of  the 
increased  TDOA  measurement  noise  are  evident  at  the  long  range  toward  the  end 
of  the  run. 

The  overall  tracking  capability  in  an  IH  medium  demonstrated  by  this  run  is 
also  shown  in  Fig.  5.6  as  am  XY  plot.  The  estimated  track  (heavy  line)  approaches 
the  actual  track  (marked  with  dots)  whereas  the  measurements  are  very  noisy  (the 
fluctuating  line). 
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D.  EFFECTS  OF  ERRONEOUS  SVP 


The  sensitivity  of  the  IH  MP  inversion  to  errors  resulting  from  an  assumed 
SVP  which  is  different  from  the  true  SVP  is  examined  in  this  section.  A  SVP  with 
positive  surface  gradient  of  0.105  sec-1  was  used  as  the  actual  medium  in  Run 
11  (Refcase  A4261).  A  similar  but  slightly  different  SVP  with  a  surface  gradient 
of  0.100  sec-1  was  used  in  the  prefilter  (Evalcase  A3253).  The  results  are  shown 
in  Fig.  5.8.  The  range  error  is  quite  small  for  ranges  less  them  7  km.  At  ranges 
greater  than  8  km  the  TDOA  resulting  from  the  true  SVP  using  the  Refcase  axe 
not  covered  by  the  direct  function  used  by  the  prefilter.  This  again  is  not  a  limit 
of  the  basic  method  but  rather  a  result  of  the  limitation  of  the  current  SMART 
model  which  forced  the  grid  selection  (see  Appendix  F).  For  the  region  evaluated 
here  one  can  conclude  that  the  method  is  insensitive  to  small  errors  in  SVP.  This 
trend  is  expected  to  extend  to  longer  ranges. 

When  large  SVP  errors  are  present,  the  resultant  tracking  error  is  also  large. 
This  is  demonstrated  by  Run  12  where  a  gradient  of  0.105  sec-1  was  used  in  the 
Refcase  and  a  gradient  of  0.05  sec-1  was  used  as  an  Evalcase  (C4261  and  C2251 
respectively).  The  results,  shown  on  Fig.  5.9,  indicate  large  tracking  errors  that 
resemble  the  effect  of  the  uncompensated  IH  medium  demonstrated  in  Run  8.  This 
is  expected  since  the  SVP  error  in  the  assumed  SVP  in  Run  12  is  of  the  same  order 
of  magnitude  as  the  toted  uncompensated  inhomogenuity  in  Run  8. 


E.  SUMMARY  OF  SIMULATION  RESULTS 


The  results  of  the  complete  system  evaluation  can  be  summarized  as  follows: 

1.  IH  propagation  introduces  significant  errors  in  MP  tracking.  The  effect  must 
be  compensated  for  in  order  to  achieve  accurate  tracking. 

2.  The  IH  prefilter  proposed  here  is  functionally  sound  and  can  provide  compen¬ 
sation  for  the  effects  of  the  IH  medium. 

3.  The  DR  grid  used  to  compute  the  direct  function  should  be  chosen  with  care. 
It  should  cover  all  possible  target  positions,  with  an  additional  margin  to 
enable  handling  of  measurement  noise.  Too  wide  a  grid,  primarily  along  the 
depth  axis,  adds  ambiguous  solutions  and  should  be  avoided. 

4.  As  a  general  rule,  the  maximum  search  a-~<*  size  should  be  made  larger  than 
twice  the  standard  deviation  of  the  depth  and  range  measurement  noise,  i.e., 
area  size  >  (2 crd)  x  (2crr). 

5.  The  maneuvering  target  tracker  can  handle  reduced  effective  sampling  rates 
resulting  from  inversion  failures.  The  original  sampling  period  was  selected 
here  as  5  sec  to  match  target  dynamics  and  to  provide  a  reasonable  TDOA 
observation  time.  Effective  sampling  rates  which  at  times  Eire  reduced  to  as 
low  as  40%  of  the  original  rate  can  still  support  reasonable  tracking. 

6.  The  inversion  is  insensitive  to  small  errors  between  the  assumed  SVP  and  the 
actual  one.  (5-10%  in  the  sound  speed  gradient) 

7.  The  overall  performance  of  the  MP  tracking  degrades  rapidly  with  increasing 
range.  This  results  from  the  vanishing  TDOA  at  longer  range  and  the  rapid 
increase  in  TDOA  noise  due  to  the  attenuated  ten-get  signal. 

8.  As  a  result  of  the  growing  noise  at  longer  ranges,  measurement  biases  develop 
which  tend  to  underestimate  the  range.  Eventually  the  MP  tracking  filter  fails 
completely  and  losses  track  of  the  target.  Reliable  tracking  can  be  expected 
up  to  ranges  of  about  15  km. 


This  chapter  concludes  the  thesis.  It  includes  a  brief  review  of  the  problem; 
a  short  description  of  the  solutions  and  a  summary  of  the  results.  It  ends  with  a 
suggestion  for  further  extension  of  the  work. 

Current  acoustic  broadband  MP  tracking  techniques,  which  have  evolved  over 
the  last  ten  years,  axe  based  on  straight  line  propagation.  In  reality  the  ocean 
medium  is  inhomogeneous  and  the  sound  propagates  along  curved  ray  paths.  MP 
tracking  therefore  encounters  large  errors  if  the  IH  propagation  is  not  accounted 
for.  Since  there  is  no  closed  form  relation  between  the  TDOA  and  source  depth  and 
range  (DR)  for  the  inhomogeneous  case,  the  compensation  for  the  inhomogeneity 
is  computationaly  complex  and  until  now  had  not  been  done. 

Previous  MP  tracking  algorithms  have  also  assumed  that  the  TDOA  are  al¬ 
ways  resolvable  and  associable  with  the  corresponding  path.  In  practical  tracking 
situations  the  TDOA  axe  not  always  resolvable  due  to  finite  bandwidth  of  the  re¬ 
ceiver  and  association  of  TDOA  with  the  corresponding  path  is  very  difficult.  The 
impact  of  these  assumptions  on  MP  tracking  is  significant. 

In  addition  TDOA  measurements  axe  noisy.  The  SNR  depends  on  the  strength 
of  the  received  signal  which  in  turn  depends  on  range  due  to  the  propagation  loss. 
A  realistic  overall  evaluation  of  a  MP  tracking  system  requires  a  more  accurate 
model  for  the  noise  than  those  previously  used,  as  well  as  a  simulation  of  the  IH 
medium  effects. 

Current  MP  tracking  algorithms  employ  the  multiple  model  (MM)  estimation 
technique.  This  technique  produces  tracking  biases  since  the  multiple  models  are 
in  most  cases  not  centered  around  a  true  model  of  the  target. 
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This  work  focused  on  developing  a  model  that  accounts  for  the  IH  propagation 
and  the  effects  of  realistic  TDOA  processing.  A  key  contribution  of  this  work  is  the 
development  of  an  algorithm  to  compute  depth  and  range  from  TDOAs  produced 
by  MP  in  an  IH  medium.  In  our  approach  the  TDOAs  are  assumed  to  be  estimated 
by  realistic  instrumentation  and  thus  have  limited  resolution  and  are  not  assumed 
to  be  associable  with  their  corresponding  path. 

The  performance  of  the  model  is  evaluated  using  a  accurate  range  dependent 
model  for  the  noise.  The  depth  range  calculation  is  performed  with  the  help  of  an 
eigenray  program  exercised  over  a  DR  grid.  The  realistic  TDOAs  are  numerically 
computed  offline  and  stored  in  a  direct  function  table.  The  measured  TDOAs 
are  then  used  online  to  produce  the  corresponding  DR  using  linear  interpolation. 
The  method  accounts  for  any  multipath  structure  and  is  not  dependent  on  the 
association  of  the  TDOA  with  the  multipath. 

The  interpolation  is  performed  as  follows.  The  DR  and  TDOA  values  of  three 
points  in  the  direct  function  table  that  have  TDOA  values  close  to  the  measured 
set,  are  used  in  the  interpolation.  The  triangle  of  the  three  closest  points  in  the 
direct  function  table  is  located  by  means  of  a  unique  search  algorithm  called  the 
local  search.  The  search  is  limited  to  a  small  area  centered  around  the  predicted 
target  position  in  order  to  improve  the  efficiency  and  reduce  the  ambiguity. 

A  technique  for  centering  the  MM  hypotheses  bank  around  the  estimated 
target  maneuver  command  to  reduce  the  tracking  bias  was  also  developed.  The 
entire  model  is  shifted  in  a  manner  that  requires  no  computation.  It  is  based  on  an 
analytically  derived  linear  relation  between  the  command  hypothesis  bank  and  the 
states  of  the  corresponding  models.  A  second  order  smoother  instead  of  the  first 
order  one  previously  used  eliminates  the  lag  error  encountered  in  earlier  tracking 
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designs.  An  evaluation  of  the  performance  of  the  new  MP  tracking  scheme  in 
shallow  waters  using  a  realistic  model  for  the  noise  shows  very  good  results. 

A  comprehensive  set  of  test  cases  were  run  to  evaluate  the  components  of 
the  MP  tracker  both  individually  and  as  a  total  system.  The  main  results  of  the 
simulation  show  that: 

1.  Significant  errors  are  introduced  in  MP  tracking  by  the  ocean  IH  acoustic 
medium. 

2.  The  new  MP  inversion  method  and  the  tracking  technique  can  compensate  for 
the  IH  medium  effects. 

3.  The  effects  of  limited  delay  resolution  and  lack  of  TD0A  and  path  associability 
are  also  accounted  for  by  the  new  method. 

4.  The  depth  measurement  is  more  affected  by  the  ambiguity  resulting  from  the 
lack  of  TDOA  path  associability  than  the  range  measurement. 

5.  Realistic  noise  limits  the  performance  of  shallow  water  MP  tracking  to  short 
ranges  (say  less  than  15  km).  The  noise  initially  introduces  estimation  biases 
and  eventually  leads  to  a  complete  loss  of  track. 

5.  At  short  ranges  the  method  is  insensitive  to  small  errors  in  SVP.  The  method 
is  more  sensitive  to  errors  in  ocean  bottom  depth.  Large  errors  in  SVP  or 
bottom  depth  introduce  large  MP  measurement  errors. 

7.  Recentering  the  MM  around  the  estimated  maneuver  command  and  the  use 
of  the  second  order  smoother  significantly  reduces  tracker  estimation  bias. 

8.  The  coordinate  decoupling  technique  currently  used  in  MP  tracking  introduces 
tracking  errors  for  nonmaneuvering  targets  around  the  CPA  at  short  ranges. 

9.  A  second  order  target  model  is  sufficient  for  maneuvering  target  tracking  if 
the  Kalman  gain  is  properly  trimmed. 

There  are  several  directions  in  which  this  work  can  be  extended.  However  the 
following  specific  direction  is  a  natural  extension  of  the  ideas  and  offers  potential 


Use  of  more  than  the  first  two  TDOAs  has  the  potential  of  reducing  the  depth 
and  range  estimation  noise.  More  importantly  use  of  additional  TDOAs  will  reduce 
the  ambiguity  in  target  position.  If,  in  addition,  information  on  the  reflections  is 
available*  and  the  signal  spectrum  is  known  then  the  complete  ACF  of  the  received 
signal  can  be  constructed  for  every  point  in  the  grid.  A  classification  method  could 
then  be  set  up  to  measure  the  closeness  of  the  measured  ACF  to  those  constructed 
for  the  points  in  the  grid.  The  conditional  probabilities  of  the  target  position  for 
any  of  the  DR  grid  points  can  be  generated  by  the  MM  tracker  and  incorporated 
into  the  classifier.  Interpolation  between  the  closer  DR  points  can  then  be  used  to 
compute  the  target  position. 

In  summary,  this  research  investigated  the  influence  of  realistic  conditions  on 
the  multipath  tracking  of  maneuvering  targets.  The  inhomogeneity  of  the  acoustic 
ocean  medium;  the  limitations  of  TDOA  estimation;  and  the  difficulty  in  associa¬ 
tion  of  TDOA  with  multipath  all  have  a  significant  effect  on  the  tracking.  A  new 
model  was  devised  to  counter  these  effects.  A  very  thorough  evaluation  of  the  model 
was  conducted  with  a  special  emphasis  on  the  reality  of  the  simulated  condition. 
The  results  indicate  that  in  general  the  MP  tracking  is  limited  to  short  ranges. 
Accurate  tracking  can  only  be  achieved  if  the  distorting  effects  of  the  medium  and 
the  delay  estimation  are  compensated  for  as  was  successfully  done  here. 


*  Such  as  when  operating  in  a  well  known  and  previously  measured  area 
of  the  ocean 


TRANSITION  MATRICES 


A  separate  model  of  the  form 


Xn  =  <f>Xn-i  +  TUn-i  +  «W'n_1 


(A.  1) 


is  used  to  describe  target  motion  along  the  X,  Y  or  range  R  axis.  The  components 
of  the  matrices  appearing  in  this  equation  are  given  below. 

*11  =  1 

*21  =  *31  =  *32  =  0 

*12  =  (1  -  e~aT)  /a 


*22  =  e' 


*13  =  [l  +  (a«/C  °T  -  ae  °'’T)  /(a  -  a*)] 


*23  =  (e~°wT  -  e~aT) 


a  -  a. 


*33  =  e 


—awT 


'SA 


Ml 


r1  =  (aT-l  +  e~aT).- 

a 

I\>  =  1  -  e-aT 


r3  =o 


*i  =  {t+  (1  -  e~aT)  -  —  •  (1  -  e-^r)l  — —  1 

{  L  a  au>  J  (a  -  £»„,  j 

=  [ar  (1  -  e--7)  -  a„  (l  -  e'"7)] 


These  parameters  are  easily  derived  from  Eq.  2.3.3  in  Ref.  6  by  the  change  of 
variables 

Un  ~  -Uk  (A.2) 

a 

Wn  «  —  Wk  (A. 3) 

aaw 

and 

aw  =  a.  (A.4) 

The  advantage  of  the  scaling  used  here  is  that  both  U  and  W  are  now  in  [m/sec],  for 
example  a  constant  process  noise  W  or  command  U  of  1  [m/sec]  in  the  X  direction 
would  each  yield  x  =  1  [m/sec]  at  the  steady  state. 

For  the  cross  range  channel  the  equation  is 

B  =  ^  •  B„~i  4-  ThUb  +  VfbWn-i  (^-5) 

where  the  components  of  the  matrices  are  expressed 
in  terms  of  the  previous  models  components  as 
^411  = 

^421  —  ^431  =  ^432  =  0 

^412  =  ^12 

^422  =  ^22 

^431  =  ^31  /  B 

^432  =  ^32 /B. 

^433  =  ^33 

r6i  =T1/R 

r62  =t2/r 
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APPENDIX  B 


CONDITIONAL  MEAN  ESTIMATE 
A  constant  U  is  to  be  estimated  from  the  observation  z  —  U  4-  n  where 
n  ~  iV(0,  o2  )•  The  optimal  mean  square  estimate  is  the  conditional  mean  U  =  E{U\z). 
Two  equiprobable  hypotheses  U\,  U2  are  assumed,  which  span  the  range  of  the  ex¬ 
pected  U.  The  estimate  is  thus  given  by 

v  =  =  ^>,-  mw  =  MUi (®-i) 

.=1  ,=1 

where  f(z)  is  the  density  function. 

With  P{U\)  =  P(U2 )  =  1/2  and  the  given  normal  distribution  of  the  noise 


the  estimate  U  is  given  by: 


(J.  •  -4 _ e— (X  — t/!  )*/2<r2  .  U  .  _  1  .  g-(i-t/j)J/2<rJ 

V2*tr  i  V2*<r 

_1  .  e— (z— l/i)J/2<rJ  +  _1  e-(z-U,y/2^ 

VZirtr  V2xir 


(B- 2) 


which  after  some  algebra  is: 


f/  = - — _ | _ _ _ 

(P.3) 

If  we  let 


U,=U2-U1;Uc  =  (U1+U2)/2; 


(BA) 


e*  =  — f(z  —  Uc) 


(P.5) 


then  we  obtain 

=  i_  U2  Uit~Ql2  (e°/2  +  e~a/2)  +  U2ta<2  (e~a/2  4-  e~Q/2) 
1  +  ea  +  1  +  e-«  (e«/*  +  e~a/2)2 


(P.6) 
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(c«/2+e-«/2) 


2  e“/2  +  e-Qr/2 


After  substitution  of  Eq.  (B.5)  the  result  is 


U  =  Uc+y  tanh  0  -  Uc) 


APPENDIX  C 
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APPENDIX  D 

DELAY  ESTIMATION  ERRORS 

The  transformation  of  the  acoustic  noise  to  a  biased  range  dependent  delay 
noise  by  the  time  delay  estimator  is  considered  here.  This  noise  modeling  was 
done  in  order  to  provide  a  more  accurate  simulation  than  previous  MP  tracking 
simulations. 


A.  VARIANCE  DEPENDENCE  ON  SNR 

The  minimum  error  variance  in  estimation  of  a  variable  from  a  measurement 
with  zero  mean  random  noise  is  given  by  the  Cramer  Rao  lower  bound  (CRLB). 
The  bound  has  been  applied  to  TDOA  estimation.  Ianniello  [Ref  14]  gives  a  cur¬ 
rent  and  relevant  review  of  time  delay  estimation  error.  The  performance  of  an 
autocorrelator  time  delay  estimator  is  evaluated  there  both  analytically  and  exper¬ 
imentally. 

Iannoello  shows  that  if  the  x(t)  signal  and  the  noise  n(t)  both  have  uniform 
spectral  density  in  the  band  0  -  B  Hz  and  are  zero  outside  of  that  range,  then  the 
variance  of  the  autocorrelation  estimate  of  r  from  the  single  multipath  signal 

y(t)  =  x(t)  +  ax(t  -  t)  -(-  n(t)  (LU) 


where  a  is  the  relative  attenuation  (gain)  of  the  path  is  given  by 

[(1  +  a2)2  4-  a2]  Sg  +  2(1  +  a2)SB  •  NB  +  N2 

—  yr?  k  - 

a2S2 

1  3 


°V  —  ao 


a0  = 


BT  v2B2 


(D.  2) 
(D.  3) 
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where  S  and  N  axe  the  spectral  densities  of  the  signal  and  the  noise  in  the 
baseband  0  —  B  and  T  is  the  observation  time.  For  the  special  case  of  a  —  1  the 
variance  reduces  to 

4  =  4  •  ■  [5  +  4  SNR->  +  SNR-2]  (DA) 

which  has  three  dominant  regions  depending  on  the  signal  to  noise  ratio 

SNR  =  S/N.  (D.5) 


For  SNR  >>  1  the  variance  approaches  a  constant  indicating  the  performance  limit 
of  this  type  of  estimator. 


= 


7r2B3r‘ 


(D.  6) 


For  SNR  <<  1  the  error  variance  is  inversely  proportional  to  the  second  squared 
SNR.  In  between  these  two  regions  there  is  a  transition  from  constant  to  dependence 
on  SNR”2. 

In  his  paper  Ianniello  indicates  close  agreement  of  Eq.  (D.2)  both  with  sim¬ 
ulation  results  and  with  the  CRLB  for  the  typical  case  of  SNR  <<  1  encountered 
in  practice. 


B.  VARIANCE  DEPENDENCE  ON  RANGE 

For  a  stationary  acoustic  noise  with  bandlimited  spectral  density  NB  and  signal 
intensity  given  by  Eq.  (3.3)  the  SNR  at  the  input  of  the  delay  estimator  becomes 
range  dependent  and  is  given  by 


SNR(R)  =  SNR0  •  Rp 


(D.  7) 


where  SNR0  is  the  SNR  at  range  R  =  1  [m]. 


Two  models  were  used  here  to  represent  the  range  dependence  of  the  noise. 
The  first*  was  a  simple  dependence  on  some  power  p  of  range  within  a  constant  of 
proportionality  <To  selected  to  match  realistic  noise  at  ranges  5  <  R  <  15  km. 

cr  =  (70  •  R?  (D. 8) 

This  model  was  used  to  develop  the  general  trend  of  the  inversion  and  the  tracker 
when  the  noise  is  range  dependent.  The  second  noise  model*  was  based  on  the 
exact  Eq.  (D.2)  with  SNR  being  the  function  of  range  given  in  Eq.  (D.6).  Some 
examples  of  the  STD  of  TDOA  as  function  of  range  for  typical  values  of  SNR  are 
shown  in  Fig.  D.l  and  D.2. 

C.  DELAY  ESTIMATION  BIAS 

Two  of  the  sources  contributing  to  MP  bias  errors  are  associated  with  the 
delay  estimation.  These  are  the  nonnegative  character  of  the  delay  estimate  and 
the  effect  of  close  ACF  peaks.  These  are  discussed  next. 

The  symmetry  of  the  autocorrelation  function,  namely 

ACF(r)  =  ACF(-r)  (D.  9) 

prevents  the  discrimination  between  positive  and  a  negative  time  delay  (a  lagging 
versus  a  leading  replica).  Since  in  the  multipath  situation,  there  really  is  a  first 
arrival  followed  by  lagging  replicas,  the  use  of  the  positive  part  of  the  time  lag  axis  is 
conceptually  justified.  This  however  implies  that  the  error  in  time  delay  estimation 
is  not  a  simple  gaussian,  but  rather  reflects  around  the  r  =  0.  A  distribution  which 
is  equal  to  0  for  negative  time  delays,  is  used  here  to  describe  the  estimate  r.  An 


*  Noise  type  (NT)  one. 

*  NT  two 


estimate  r.  An  estimate  with  this  distribution  is  generated  by  taking  the  absolute 
value  of  the  sum  of  the  time  delay  and  a  normal  zero  mean  noise. 

When  the  actual  time  delay  is  large  relative  to  the  estimation  error  variance, 
the  effect  of  reflection  around  the  0  lag  is  marginal.  When  the  actual  delay  is  small 
compared  to  the  error  variance  the  resulting  error  bias  has  a  stronger  effect. 
Figure  D.3  shows  three  distributions  of  the  process 

t  =  |<  +  n|  (Z?.10) 

for  t  taking  the  value  of  1.8<7, 0.85(7, 0.22(7  w^er"  n  ~  jV( 0,  a2).  It  is  evident  that 
when  t  is  smaller  than  0.5<r  it  will  be  significantly  overestimated  which  will  cause 
an  underestimated  range. 


Fig.  D.3.  Nonnegative  bias. 


D.  CLOSE  ACF  PEAKS 


When  the  actual  delays  are  short  compared  to  the  signal  correlation  time  ( 1  /B) 
the  presence  of  an  adjacent  peak  will  shift  the  location  of  the  estimated  peak.  This 
shift  can  be  either  positive  or  negative  depending  on  the  specific  set  of  MP  and 
resulting  ACF.  Since  it  is  analytically  predictable  it  can  conceptually  be  corrected 
for  as  is  done  in  the  maximum  likelihood  estimator  [14].  This  error  was  therefore 
not  simulated  here. 


APPENDIX  E 


DETAILED  DATA  DEFINITION  AND  SELECTION 

This  appendix  describes  in  detail  the  parameters  used  in  the  various  simula¬ 
tions  and  evaluations.  Run  data  parameters  axe  described  in  Section  A.  The  values 
selected  for  the  10  simulation  Runs  are  listed  in  Section  B.  Section  C  describes  the 
case  code  name. 

A.  SIMULATION  RUN  DATA 
Target  and  Scenario  Data 

Ad  -  [sec]  The  target’s  depth  time  constant  (1  /aj,  in  Eq.  (2.5). 

Ar  -  [sec]  The  target’s  horizontal  time  constant  (1  /a  in  Eq.  (2.2). 

Pin  -  The  pilot  number  1,  2  or  3  per  Chapter  Two,  Section  E. 

Medium  and  Prefilter 

Case  -  A  selection  of  parameters  defining  the  ocean  medium 
and  its  acoustic  modeling. 

Mdn  -  Medium  number  1  through  4. 

Pfn  -  Prefilter  number  1  through  6. 

p  -  [number]  Power  of  the  range  dependence  of  the 

delay  estimation  variance  (See  Chapter  Two, 

Section  D.4,  Chapter  Three,  Section  A,  and  Appendix  D.) 

P6  -  [number]  Power  of  range  dependence  of  bearing 
estimation  noise,  (and  of  DR  noise  when  it  is 
used  with  medium  and  prefilter  number  1,  Eq.  (2.101)). 
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Two  cases  were  used  for  each  Run. 

Refcase  -  The  reference  is  a  set  of  parameters  considered  as 
the  “actual”  conditions  used  by  the  medium 
program  to  convert  the  actual  DR  position  to  TDOA. 

Evalcase  -  The  evaluated  ocean  conditions  and  used  by 
the  prefilter  in  the  inversion  of  TDOA  to  DR. 

NT  -  Noise  type  1  =  the  simplified  model  Eq.  (2.96). 

2  =  the  improved  model,  Appendix  D. 

crtl  -  [ms]  Standard  deviation  of  TDOA 

t\  (or  T\ )  estimate  at  range  of  500  m.  This  parameter 

is  was  used  only  by  noise  model 

type  1,  but  is  computed  for  noise  model 

type  2  for  reference  as  well. 

<7(2  -  Same  as  above  but  for  <2 (r2)- 

<7 d  -  [m]  Depth  error  standard  deviation  at  range  of 

500  m,  used  with  medium  1  and  prefilter  1. 

<7 r  -  [m]  Range  error  standard  deviation  at  range  of 

500  m,  used  with  medium  1  and  prefilter  1. 

<7^  -  [deg]  Bearing  standard  deviation  at  range  of  500  m. 

Bw  -  [Hz]  Receiver  bandwidth. 

Atn  -  [dB]  A  vector  of  surface  and  bottom  bounce  losses. 

SNRo  -  [dB]  Signal  to  noise  ratio  for  target  at  range  of  1  m. 

measured  at  the  input  to  the  delay  estimator. 

MM  Parameters 

T  -  [sec]  System  sampling  time  interval. 

Ord  -  MM  target  model  order,  2  or  3. 

Aq  -  [sec]  Command  noise  coloring  time  constant  (1  /aw  Eq.  (2.7)). 
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Umin  -  [m/sec]  Lowest  speed  command  in  the  hypothesis  bank. 

Umaz  -  [m/sec]  Highest  speed  command  in  the  hypothesis  bank. 

N  -  Number  of  filters  in  the  MM. 

Pr  -  Command  transition  probability  0  to  1. 

Anl  -  W  averaging  constant,  (iterations  to  reach  63%)  Eq.  (2.23b) 

An3  -  W 0p,Uop  averaging  constant,  (iterations  to  reach  63%) 

Eq.  (2.26)  and  (2.27) 

An4  -  Xop  averaging  constant,  (iterations  to  reach  63%) 

Eq.  (2.28) 

An  -  Uc  averaging  constant,  (iterations  to  reach  63%) 

Eq.  (2.88) 

B.  Runs 

The  specific  values  of  the  parameters  for  each  run  are  given  below: 

Run  -  1 

Target  and  scenario  data 

Ad  -  35  sec. 

Ar  -  40  sec. 

Pin  -  3 

Mdn  -  1 

Pfn  -  1 

Pb  -  1 

(J  d  -  1  [m] 

(Jr  -  50  [m] 

(Jb  -  1  [deg] 


MM  parameters 


Aw 

-  25  [sec] 

Umin 

-  —30  [m/sec] 

U max 

-  +30  [m/sec] 

N 

-  11 

Pr 

-  0.9 

An  1 

-  10 

An3 

-  4 

An4 

-  5 

An9 

-  100 

Run 


Pin 

Mdn 

Pfci 

Refcase 

Evalcase 

NT 

P 

Pb 

On 

at2 

Ob 

Bw 

Atn 

SNR 

Run 


P 

on 

Ot2 

SNR 

Run 


P 


2  (All  parameters  are  the  same  as  for  Run  1 
except  the  following:) 

3 

3 

3 

sun 

sim 

2 

1 

1 

0.005  [m] 

0.004  [m] 

0.08  [deg] 

2000  [Hz] 

5.5  [dB] 

50  [dB] 

3  (all  the  parameters  are  the  same  as  for  Run  2 
except  the  following:) 

2 

0.05  [ms] 

0.05  [ms] 

70  [dB] 

4  (All  parameters  are  the  same  as  Run  2  except 
for  the  following:) 

2 


US? 


<7(2 

SNR 

Run 


P 

<7*1 

<7(2 

Atn 

Run 


Pin 


-  0.08  [ms] 

-  60  [dB] 

-  5  (All  parameters  Eire  the  same  as  for  Rim  2 
except  the  following:) 

-  2 

-  0.018  [ms] 

-  0.03  [ms] 

-  5.10  [dB] 

-  6  (All  parameters  are  the  same  as  Run  2 
except  for  the  following:) 

-  1 


Mdn 

Pfh 

P 

Pb 

<rti 

<7(2 

max  area  size 

Atn 

SNR 


4 

5 
2 
1 

0.018  [ms] 
0.018  [ms] 
8x15 
5.5  [dB] 

50  [dB] 


Run 


7  (All  parameters  are  the  same  as  Run  6 
except  for  the  following:) 


max  area  size  -  5x9 


Run 


Refcase 

Evalcase 

NT 

P 

Pb 

<7(1 

<7(2 


8  (All  parameters  are  the  same  as  Run  6 
except  for  the  following:) 

C2251 

Sllll 

2 

2 

2 

0.005  [ms] 

0.005  [ms] 
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55  [dB] 
70  [dB] 
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Atn 

SNR 


Run 


9  (All  parameters  are  the  same  as  in  Run  8 
except  for  the  following:) 


Atn 

SNR 


5.10  [dB] 
50  [dB] 


Run 


10  (All  parameters  axe  the  same  as  in  Run  9 
except  for  the  following:) 


Refcase 

Evalcase 


C2251 

C2251 


Run 


11  (All  parameters  are  the  same  as  in  Rim  10 
except  for  the  following:) 


Refcase 

Evalcase 


A426 

A3253 


Run 


12  (All  parameters  are  the  same  as  in  Run  11 
except  for  the  following:) 


Refcase 
Evalcase 
<*ti 
at  2 

Atn 

SNR 


C4261 
C2251 
0.005  [ms] 
0.005  [ms] 
5.5  [dB] 

70  [dB] 


C.  Case  Definition 


A  case  code  is  comprised  of  one  character  and  four  digits,  (for  example  C2461 ). 
It  was  used  to  define  the  basic  MP  assumptions.  These  include  homogeneous  versus 


IH  propagation,  the  ocean  conditions,  the  estimation  resolution,  The  DR  grid  data 


and  some  additional  parameters.  The  code  is  used  in  the  computation  of  TDOA 


surfaces  over  the  depth  range  (DR)  grid. 


R 


The  parameters  are  defined  as  follows.  The  leading  character  indicates  the 


type  of  case. 

U  -  Homogeneous  propagation,  perfectly  resolvable  and  identifiable  MP.  This 
is  an  idealized  case  where  Eq.  (1.2)  is  used  to  compute  the  unsorted  (hence  the  U) 
time  delays  tit%.  It  is  used  mainly  for  methodology  purposes,  demonstrating  the 
differences  between  the  TDOA  and  modified  TDOA. 

S  -  Homogeneous  propagation,  perfectly  resolvable  but  non  identifiable  paths. 
Eq.  (1.2)  is  again  used  to  generate  the  time  delays  which  are  then  sorted  according 
to  delay  magnitude  (hence  the  S)  at  each  DR  point  and  denoted  by  tj  and  *2- 
The  variable  t\  is  now  always  the  shortest,  but  it  cam  no  longer  be  associated  with 
surface  reflection  nor  can  the  longer  t<i  be  associated  with  the  bottom  bounce.  This 
case  is  useful  as  an  idealized  reference  for  the  later  IH  limited  resolution  since  they 
are  conceptually  similar. 

A,B,C,D  -  The  realistic  IH  limited  resolution  non  identifiable  cases.  The  letter 
indicates  the  receiver  time  delay  resolution  in  milliseconds  as  follows. 

A  -  0.05  C  -  0.5 

B  -  0.2  D  -  1.0 

The  four  digits  define  the  parameters  required  for  the  numerical  generation  of 
the  TDOA  surfaces  as  follows.  The  first  digit  indicates  the  SVP.  Four  SVP’s  of 
varying  inhomogenuity  are  used  throughout  this  research  as  shown  below: 
case  1  st  digit  1.  2.  3.  4. 


■  » 'w  j.>  .  ,>fc»  .  *  .  ‘•t  'j*  "  »  '  , 


The  SVP  digit  is  of  course  meaningless  for  the  S  and  U  cases  where  homage 
neous  propagation  (SVP  digit  =  1)  is  always  assumed. 

The  second  digit  details  the  observer  and  ocean  depth  information. 


TABLE  E.l 

DEPTH  DATA 

case  2nd  digit 

1 

2  3 

4 

5 

Ocean  depth  m 

503 

513 

518 

550 

Observer  depth 

162 

162 

162 

162 

The  third  and  fourth  digits  describe  the  discrete  depth  and  range  grid  used 
for  the  direct  function  table. 

TABLE  E.2 
GRIP  DATA 


Grid  density 

1 

2 

3 

4 

5 

Initial  depth  [m] 

160 

- 

160 

160 

160 

Depth  Step  size  [m] 

-50 

25 

-100 

10 

Number  of  steps 

11 

21 

6 

11 

Initial  range  [m] 

500 

500 

500 

500 

500 

Range  step  size  [m] 

500 

1000 

500 

1000 

250 

number  of  steps 


25 


25 


49 


13 
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APPENDIX  F 


THE  EIGENRAY  ACOUSTIC  MODEL 

This  appendix  describes  the  selection  of  the  eigenray  model  and  the  evaluation 
of  the  SMART  model. 

A.  REQUIREMENTS  AND  MODEL  SELECTION 

The  general  goal  of  the  eigenray  model  is  to  find  and  trace  all  the  rays  passing 
‘1 .  ^ugh  the  two  specified  receiver  and  source  end  points.  The  following  is  required 
from  an  eigenray  model  if  it  is  to  be  used  for  MP  tracking. 

1.  Correctness  and  High  Accuracy 

The  total  sound  travel  time  from  the  source  to  the  receiver  may  be  of 
the  order  of  a  few  tens  of  seconds.  The  milliseconds  of  travel  time  difference, 
while  vital  for  the  MP  ranging,  are  a  very  small  percent  of  the  total  travel  time. 
This  is  especially  true  for  the  longer  ranges  where  the  travel  time  grows  while  the 
differences  decrease.  The  resulting  output  accuracy  requirement  is  of  the  order  of 
.001%  which  dictates  even  higher  precision  in  the  model’s  internal  computations. 
The  issue  of  correctness  relates  to  the  possibilities  of  not  finding  an  existing  path 
or  producing  imaginary  nonexisting  paths;  it  too  is  related  to  the  computer  word 
length. 

2.  Compactness 

The  eventual  use  of  a  MP  tracking  algorithm  in  a  real  time  systems  appli¬ 
cation  imposes  additional  constraints  on  run  time,  supporting  hardware,  and  size 
of  the  model.  Research  oriented  programs  requiring  extra  long  computer  words  or 
computation  time  will  not  fit  the  smaller  data  processors  usually  available  for 


target  tracking  in  onboard  systems.  What  is  required  is  a  program  running  on  a 
microcomputer  that  can  solve  a  typical  eigenray  problem  in  a  matter  of  seconds. 


B.  MODEL  SELECTION 

Some  eigenray  models  were  developed  in  recent  years,  but  the  conflicting  re¬ 
quirements  of  accuracy  and  compactness,  narrows  the  selection  to  two  programs. 
One  model  called  SMART,  for  SMall  Acoustic  Ray  Tracer  was  developed  by  Novick 
of  Mission  Sciences  Corp.  [Ref.  22];  the  other  is  a  research  tool  developed  by  Spies- 
berger  as  an  application  of  an  earlier  model  [Ref.  31]. 

The  SMART  model  was  selected  because  it  was  better  integrated  as  a  working 
model  and  provided  more  detailed  output  data.  The  Spiesberger  model  has  is  the 
advantage  that  it  can  account  for  a  sloped  ocean  bottom.  However  this  model  lacks 
a  reliable  transmission  loss  output  and  is  not  well  integrated  for  the  intended  use 
here. 

The  eigenray  model  is  used  in  a  fashion  that  is  separable  from  the  rest  of  the 
MP  inversion  or  target  tracking  algorithm.  This  provides  for  easy  replacement  or 
upgrading  of  the  model.  This  also  led  to  our  decision  to  use  the  more  reliable 
integrated  model  for  this  first  IH  compensation  attempt,  and  maintain  the  option 
to  upgrade  and  extend  it  in  the  future. 

The  SMART  model  was  evaluated  extensively  and  the  main  result  of  this 
evaluation  is  presented  in  Section  D. 

C.  SMART  MODEL  -  A  BRIEF  DESCRIPTION 

Finding  the  eigenrays  is  the  first  task  accomplished  by  the  model.  Source  and 
receiver  depths  are  set  to  the  specified  input  parameters  and  an  iterative  search 


for  the  initial  angle  of  the  ray  which  will  reach  the  receiver  depth  at  the  specified 
range,  is  conducted.  The  interval  of  possible  ray  angles  is  divided  into  sectors  for 
which  the  dependence  of  range  on  the  angle  is  a  continuous  well-behaved  monotonic 
function. 

Eigenrays  are  searched  for  in  every  sector  that  spans  the  desired  range  by 
successive  approximation,  the  procedure  terminates  when  the  eigenray  reaches  a 
specified  range  accuracy.  Only  rays  with  a  total  number  of  bounces  less  them  a 
specified  input  value  are  considered.  Forward  ray  tracing  equations  similar  to  those 
presented  in  Chapter  Three  are  used  to  locate  the  limits  of  the  angle  sector,  to  try 
out  angles  in  the  search  process,  and  to  compute  the  detailed  output  for  every 
eigenray  found. 

D.  LIMITATION  OF  THE  MODEL 

The  program  does  not  find  eigenrays  that  are  within  ±0.15  deg  from  the 
horizontal.  This  is  partially  caused  by  the  limitation  of  intensity  computation  at 
these  angles  (recall  discussion  in  Chapter  3  and  Eq.  (3.20)).  Another  cause  for 
the  limit  is  the  finite  computer  word  length  and  the  increment  used  for  searching 
the  ray  around  the  horizontal  direction.  The  ±0.15°  limitation  relates  to  both 
angle  at  the  source  and  at  the  receiver.  An  example  of  the  problem  is  shown  in 
Fig  F.l.  The  bounce  count  (BC)  of  case  Cllll*  computed  using  SMART  shows 
a  high  value  along  the  receiver  depth  at  long  range.  No  loss  of  a  path  or  loss  of 
resolution  should  give  rise  to  this  phenomena.**  The  unexpected  BC  was  traced 
*  See  Appendix  E  for  case  definition. 

**  Near  the  surface  and  the  bottom  the  corresponding  single  bounce  is  not 
resolvable  from  the  direct  path.  At  depth  of  —170  m  the  single  surface  and 


of  the  eigenray  model  to  detect  the  direct  path.  This  results  from  the  direct  ray 
being  close  to  horizontal  at  those  ranges.  (The  vertical  separation  between  the 
grid  depth  line  at  152  m  and  receiver  at  162  m  is  10  m.)  These  limitations  were 
brought  to  the  attention  of  the  author  of  the  SMART  model  and  axe  currently 
being  corrected.  Since  the  model  is  employed  in  the  prefilter  in  a  modular  way, 
these  limitations  do  not  effect  the  generality  of  the  proposed  inversion  procedure. 

The  choice  of  The  SVP  and  geometry  selected  for  use  in  this  work  was  de¬ 
signed  to  avoid  this  problem  by  limiting  the  grid  to  a  region  were  the  rays  are  not 
horizontal  at  all.  This  was  done  by  (1)  restricting  the  DR  grid  to  +60  -  +160  m 
while  the  receiver  was  placed  at  162  m;  (2)  limit  the  maximum  range  to  25  km; 
and  (3)  restricts  the  SVP  gradients  to  positive  surface  gradients  only.  The  overall 
geometry  is  shown  in  Fig.  F.2. 


bottom  bounce  are  not  resolvable. 
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APPENDIX  G 

THE  LOCAL  SEARCH  ALGORITHM 


In  this  Appendix  the  component  of  the  local  search  algorithm,  the  inversion 
test,  the  inward  search  and  the  outward  serrch,  are  described.  A  pseudo  code 
description  of  the  algorithm  is  also  included.  The  input  and  output  parameters 
used  by  the  algorithm  are  defined  in  Chapter  Three. 

A  Numerical  analysis  procedure,  developed  for  finding  zeros  of  complex  func¬ 
tions  [Ref.  29]  was  modified  to  apply  to  the  search.  For  this  purpose  the  search  has 
to  be  viewed  in  the  context  of  the  graphical  statement  of  the  problem  as  discussed 
in  Chapter  Three,  that  is,  as  finding  the  DR  coordinates  of  the  point  of  intersection 
of  the  constant  and  <2  contours  (see  Fig.  3.12).  The  inversion  problem 


can  be  rewritten  as 


(G.l) 


Where  the  D  and  R  are  the  unknowns.  If  the  vector  function  /  is  written  as  two 


separate  functions  f\  and  /2  Eq.  (G.2)can  be  written  as  the  two  equations 


d-^C;)=° 


(G.3a) 


(G.3b) 


The  contour  lines  of  Fig.  3.12  represent  loci  of  the  DR  solution  of  Eq.  (G.3). 
Every  point  p  on  the  DR  grid  with  depth  Dp  and  range  Rp  can  be  defined  in  terms 
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of  it’s  relation  to  the  input  t\  contour  by  a  one  digit  number  Cpi  given  by 


f  l  ;tip  >  1 1 

Cpl  =  <  0  ;  tlp  =  *,  (G.4) 

l  1  ;  tip  <  ti 

and  in  terms  of  the  it’s  relation  to  the  second  input  contour  <2  by  another  digit 
CP2  given  by 

f  1  \hp>t2 

Cp2  =  {  0  ;f2p  =  22  (G.5) 

l  —  1  ;*2 p  <  <2 

and  CP2  are  referred  to  as  characteristic  digits.  An  example  of  the  two  char¬ 
acteristic  digits  for  some  points  is  shown  on  Fig.  G.l. 

A.  THE  INVERSION  TEST 

A  1  x  1  rectangle  can  contain  a  solution  only  if  its  characteristic  digits  c,i 
and  CP2  are  either  zero  or  if  they  acquire  at  least  two  different  values  in  one  of  the 
four  rectangle  corners. 

Some  examples  are  shown  in  Fig.  G.l(c).  In  (a)  both  Cp\  and  Cp2  are  1  and 
—  1  at  different  points.  In  (b)  Cp\  is  0  and  1,  while  C? 2  is  1  and  —1.  This  rule  is  a 
simplified  version  of  a  rule  analytically  derived  in  Ref  [29]  from  Cauchy’s  principal 
value  argument. 

Note  that  this  simple  test  is  a  necessary  but  not  a  sufficient  condition.  Fig.  G.lc 
shows  an  example  of  a  rectangle  which  succeeds  the  test  but  does  not  contain  a 
solution.  If  this  simple  conditioned  test  succeeds  it  is  followed  by  a  second  more 
rigorous  test  which  checks  to  see  if  the  two  contours  actually  intersect  or  just  pass 
through  the  rectangle.  This  is  done  by  actually  attempting  to  locate  the  solution, 
by  means  of  the  inward  seaxch  described  later. 


B.  THE  OUTWARD  SEARCH 


The  outward  search  is  a  recursive  procedure.  It  performs  the  inward  search 
over  successively  growing  search  rectangles.  If  a  call  to  the  inward  search  returns 
in  a  failure  (no  solutions  in  the  current  rectangle)  the  search  area  is  increased  by 
tripling  the  sides  of  the  current  rectangle.  A  new  current  rectangle  is  thus  formed 
centered  around  the  previous  one  with  an  area  nine  times  larger.  The  outward 
search  is  now  called  recursively  to  search  the  new  rectangle.  Again  it  tests  for 
existence  of  a  solution  by  calling  the  inward  search.  Note  that  in  its  growth  the 
outward  search  may  in  one  step  include  two  or  more  solutions  which  would  then 
all  be  located  by  the  inward  search. 

The  search  rectangle  area  is  increased  by  a  9.T  ratio  by  means  of  tripling  the 
length  of  both  its  sides.  This  is  done  in  order  to  minimize  the  duplicated  effort 
associated  with  re-  examining  the  original  rectangle  as  part  of  the  recursive  search 
of  the  new  rectangle.  However  the  rectangle  is  increased  only  up  to  the  portion 
of  the  maximum  search  area  size  that  still  remains  within  the  limits  of  the  grid. 
If  the  inward  search  for  the  current  rectangle  is  unsuccessful  and  the  search  area 
cannot  be  further  increased  in  any  direction,  the  procedure  terminates  in  a  failure, 
without  any  solution. 

Upon  termination  the  search  procedure  outputs  the  size  of  the  current  area, 
the  number  of  solutions  found  and  the  solution  quad  triangles.  If  more  than  one 
solution  is  found  the  calling  inversion  procedure  will  select  the  most  likely  one,  i.e., 
the  one  closer  to  the  predicted  position  after  the  interpolation  is  performed.  An 
example  of  the  growing  search  rectangles  is  shown  on  Fig.  3.18  in  the  dashed  line. 
The  search  is  initiated  at  point  Dv ,  f?p,  and  the  existence  of  solution  is  detected 
after  two  iterations. 


C.  THE  INWARD  SEARCH 

The  inward  search  attempts  to  locate  a  solution  in  a  given  input  rectangular 
search  area.  It  starts  by  performing  the  conditional  inversion  test  on  the  complete 
area.  If  the  test  succeeds,  the  region  is  divided  into  four  (quadrant)  rectangles 
by  halving  both  of  its  axes.  Each  quadrant  is  then  searched  for  solutions  by  a 
recursive  call  to  the  inward  search*.  If  the  possible  existence  of  a  solution  in  the 
quadrant  is  indicated  by  the  inversion  test  the  quadrant  is  subdivided  and  the 
process  continues  recursively.  If  the  test  of  the  quadrant  fails,  the  search  in  it 
terminates.  The  inward  recursion  thus  terminates  when  either  the  search  area 
reaches  the  minimum  rectangle  (of  size  1  x  1),  or  if  no  solution  exists  in  the  a 
specific  branch  (of  the  B  tree  with  branching  factor  4). 

The  triangle(s)  in  the  lxl  rectangle  is  (are)  found  by  dividing  the  rectangle 
into  it’s  two  DR  grid  triangles  and  performing  a  test  on  each  to  see  if  it  surrounds 
the  point.  The  surrounding  test  uses  vector  product  to  verify  that  the  input  ti,t2 
point  is  on  the  same  side  of  all  lines  that  make  the  sides  of  the  triangle.  The  test  is 
as  follows:  Three  vectors  Si,  S2  and  S3  are  defined  in  the  <1  #2  plane  as  the  three 
sides  of  the  triangle  in  a  cyclic  order  (see  Fig.  G.2).  Three  corresponding  vectors 
V\,  V2  and  V3  axe  formed  by  the  input  t\  and  t2  point  and  the  three  angles  of  the 
triangle  Ai,A2,A3  respectively.  The  side  of  the  input  point  relative  to  a  “side 
vector”  S  is  given  by  the  sign  of  the  vector  product  of  the  corresponding  vectors 
5  and  V  vectors.  An  input  point  inside  a  triangle  will  be  on  the  same  side  of  all 
the  three  side  vectors.  That  is 

sign  [Vj  x  Si]  =  sign  [V2  *  S2]  =  [^3  •  S3]  (G.6) 

*  The  search  is  therefore  a  B  tree  search  with  a  branching  factor  of  4. 
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Fig.  G.2.  Vectors  for  inside  test. 


In  Fig.  G.2  the  six  vectors  axe  plotted  for  a  point  inside  (pi  and  full  line)  and  for 
a  point  outside  (pi  and  dashed  lines)  a  triangle.  The  3igns  marked  along  the  V 
vectors  axe  the  signs  of  the  corresponding  vector  product  (Vi  x  5,). 

The  solutions  from  all  the  branches  of  the  recursion  are  reported  by  the  inward 
search.  This  ensures  that  all  solutions  found  in  the  searched  axea  centered  around 
the  reference  point  are  reported  for  further  consideration.  In  the  example  shown 
in  Fig.  3.18  three  iterations  of  the  inwaxd  search  are  required  to  locate  the  actual 
triangles  that  surround  the  point.  In  this  case  there  were  two  solutions. 

Pseudo  code  of  the  local  search  algorithm  follows. 


D.  PSEUDO  CODE  FOR  THE  LS  A 

TYPE  box  :  a  square  area  in  the  depth  range  grid  defined  by 
the  DR  indexes. 

FUNCTION  inside  (  tl,  t2  ,tri) :  boolean; 

inside  1  if  the  point  tl,t2  is  inside  the  triangle  defined 
by  the  three  quads  of  tri  on  the  tlt2  plane; 

0  otherewize. 

FUNCTION  close  (  INVAR:  depth,  range) :  box 

close  the  i  j  (depth  range  )  indexes  of  the  qu.  u  box  on 
the  depth  range  grid  surrounding  the  input 
depth  range  point. 

ENDFUNCTION ; 

PROCEDURE  inversion_test  (  INVAR  :  tl,  t2 ,  current; 

OUTVAR  :  result  ) 

IF 

(CPI  equals  0  or  changes  sign  in  the  current  area) 
AND  (CP2  equals  0  or  changes  sign  in  the  current  area) 
THEN  result  :-l  ELSE  result  0; 

PROCEDURE  owtward  (  INVAR  :  tl,  t2 ,  current  ,  maxsize; 

OUTVAR  :Qlist,  number_of_solut ions , areasize) 

BEGIN 

IF  current  =  previous  THEN  exit  .’reached  the  size  limit. 

ELSE 

BEGIN 

inward  (  tl,  t2 , current , qlist,  num_of_solutions) ; 
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IF  number_of_solutions  =/=  0  THEN  exit ; Solution  found. 
Else 

BEGIN  ;  Recourse, 

current  triplsized  (current), 
outward  ( tl , t2 ,  current ,  nuxn_of_solution ,  areasize) 
END; 

ENDIF; 

END; 

ENDIF; 

ENDFROCEDURE . 

PROCEDURE  inward  (INVAR  :  tl, t2 , current 

OUTVAR  :  qlist,num_of_solutions) ; 

BEGIN 

IF  current  =1x1  THEN 

split  current  box  into  triangles  tril,  tri2; 
FOR  tri  IN  (tril,  tri2)  DO 

IF  inside  (tl,t2,tri)  THEN 
add  tri  to  qlist; 

num_of_solutions : -  num_of_solutions+l ; 
ENDIF; 

ENDFOR; 

ELSE 

qlist  empty; 
num_of_solution  0; 

IF  not(inversion_test (tl,  t2,  current))  THEN  exit 
ELSE 
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split  current  into  four  quadrants  1  thru  4; 
FOR  next  IN  quadrants  1  through  4  DO 


inward  (  tl, t2, next, list, num) ; 
qlist  qlist  +  list; 

num_of_solution  num_of_solution  +  1; 
ENDFOR; 

ENDIF; 

ENDIF ; 

ENDPROCEDURE. 


PROCEDURE  inversion  (  INVAR  : tl ,  t2,  DRref ,maxsize 

OUTVAR  :  depth,  range, aver_num, aver_areasize) 
var  current  :  box; 

solutionlist  :  list  of  depth  range  solutions; 

BEGIN 

current  close (  DR  refpoint) ; 

outward  (  tl,t2, current, qlist, num_of_solutions,areasize) ; 

IF  num_of_solution  =  0  THEN  depth,  range  DRref  ;a  failure 
ELSE 
BEGIN 

aver_nuxn  average_nun  x  a  +  num_of_solutions  x  (1-a)  ; 
aver_areasize  aver_areasize  x  a  +  areasize  x  (1  -  a) ; 

FOR  tri  IN  qlist  DO 

interpolate  depth  range  using  tri; 
append  solution  to  solution  list; 

ENDFOR; 

ENDIF; 

select  depth  range  solution  closest  to  DRref  point  as  depth 
range  values  to  output. 


ENDPROCEDURE . 
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